{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "MCMC.ipynb",
      "version": "0.3.2",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Q4qew2DscIPu",
        "colab_type": "text"
      },
      "source": [
        "# 马尔可夫链蒙特卡洛方法"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qKLGsUc-cNIM",
        "colab_type": "text"
      },
      "source": [
        "蒙特卡洛法（Monte Carlo method) ， 也称为统计模拟方法 (statistical simulation method) ， 是通过从概率模型的随机抽样进行近似数值计\n",
        "\n",
        "算的方法。 马尔可夫链陟特卡罗法 (Markov Chain Monte Carlo, MCMC)， 则是以马尔可夫链 (Markov chain）为概率模型的蒙特卡洛法。\n",
        "\n",
        "马尔可夫链蒙特卡罗法构建一个马尔可夫链，使其平稳分布就是要进行抽样的分布， 首先基于该马尔可夫链进行随机游走， 产生样本的序列，\n",
        "\n",
        "之后使用该平稳分布的样本进行近似数值计算。\n",
        "\n",
        "Metropolis-Hastings算法是最基本的马尔可夫链蒙特卡罗法，Metropolis等人在 1953年提出原始的算法，Hastings在1970年对之加以推广，\n",
        "\n",
        "形成了现在的形式。吉布斯抽样(Gibbs sampling)是更简单、使用更广泛的马尔可夫链蒙特卡罗法，1984 年由S. Geman和D. Geman提出。\n",
        "\n",
        "马尔可夫链蒙特卡罗法被应用于概率分布的估计、定积分的近似计算、最优化问题的近似求解等问题，特别是被应用于统计学习中概率模型的学习\n",
        "\n",
        "与推理，是重要的统计学习计算方法。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jE5JCrmed0C4",
        "colab_type": "text"
      },
      "source": [
        "一般的蒙特卡罗法有**直接抽样法**、** 接受-拒绝抽样法**、 **重要性抽样法**等。\n",
        "\n",
        "接受-拒绝抽样法、重要性抽样法适合于概率密度函数复杂 （如密度函数含有多个变量，各变量相互不独立，密度函数形式复杂），不能直接抽样的情况。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2ZBr1H0Md0gv",
        "colab_type": "text"
      },
      "source": [
        "### 19.1.2 数学期望估计\n",
        "\n",
        "一舣的蒙特卡罗法， 如直接抽样法、接受·拒绝抽样法、重要性抽样法， 也可以用于数学期望估计 （estimation Of mathematical expectation)。\n",
        "\n",
        "假设有随机变量$x$， 取值 $x\\in X$, 其概率密度函数为 $p(x)$, $f(x)$ 为定义在 $X$ 上的函数， 目标是求函数 $f(x)$ 关于密度函数 $p(x)$ 的数学期望 $E_{p(x)}[f(x)]$。\n",
        "\n",
        "\n",
        "针对这个问题，蒙特卡罗法按照概率分布 $p(x)$ 独立地抽取 $n$ 个样本$x_{1}, x_{2},...,x_{n}$，比如用以上的抽样方法，之后计算函\n",
        "\n",
        "数$f(x)$的样本均值$\\hat f_{n}$  \n",
        "\n",
        "$\\hat f_{n} = \\frac{1} {n}\\sum_{i=1}^{n}f(x_{i})$\n",
        "\n",
        "\n",
        "作为数学期望$E_{p(x)}[f(x)]$近似值。\n",
        "\n",
        "根据大数定律可知， 当样本容量增大时， 样本均值以概率1收敛于数学期望：  \n",
        "\n",
        "$\\hat f_{n} \\rightarrow E_{p(x)}[f(x)], n \\rightarrow \\infty $\n",
        "\n",
        "这样就得到了数学期望的近似计算方法：  \n",
        "\n",
        "$E_{p(x)}[f(x)] \\approx \\frac{1} {n}\\sum_{i=1}^{n}f(x_{i})$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZtV3LYrUh9jG",
        "colab_type": "text"
      },
      "source": [
        "### 马尔可夫链 \n",
        "\n",
        "考虑一个随机变量的序列 $X = {X_{0}, X_{1},..., X(t),...}$ 这里 $X_{t}$，表示时刻 $t$ 的随机变量， $t = 0, 1, 2...$. \n",
        "\n",
        "每个随机变量 $X_{t}(t=0,1,2,...)$ 的取值集合相同， 称为状态空间， 表示为$S$.  随机变量可以是离散的， 也可以是连续的。\n",
        "\n",
        "以上随机变量的序列构成随机过程（stochastic process)。\n",
        "\n",
        "假设在时刻 $0$ 的随机变量 $X_{0}$ 遵循概率分布 $P(X_{0}) = \\pi$，称为初始状态分布。在某个时刻 $t>=1$ 的随机变量 $X_{t}$与前\n",
        "\n",
        "一个时刻的随机变量 $X_{t-1}$ 之间有条件分布 $P(X_{t}|X_{t-1})$ 如果 $X_{t}$ 只依赖于 $X_{t-1}$, 而不依赖于过去的随机变量 \n",
        "\n",
        "${X_{0}，X_{1},...，X_{t-2}}$ 这一性质称为马尔可夫性，即  \n",
        "\n",
        "$P(X_{t}|X_{0},X_{1},...,X_{t-1}) = P(X_{t}|X_{t-1}), t=1,2,...$\n",
        "\n",
        "具有马尔可夫性的随机序列$X = {X_{0}, X_{1},..., X(t),...}$称为马尔可夫链， 或马尔可夫过程（Markov process)。 条件概率分布 \n",
        "\n",
        "$P(X_{t}|X_{t-1})$ 称为**马尔可夫链的转移概率分布**。 **转移概率分布决定了马尔可夫裢的特性**。\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BL8oLbgBttgs",
        "colab_type": "text"
      },
      "source": [
        "#### 平稳分布  \n",
        "\n",
        "设有马尔可夫链$X = {X_{0}, X_{1},..., X(t),...}$，其状态空间为 $S$,转移概率矩阵为 $P=(p_{ij})$， 如果存在状态空间 $S$ 上的一个分布  \n",
        "\n",
        "$\\pi = \\begin{bmatrix}\n",
        "\\pi_{1}\\\\ \n",
        "\\pi_{2}\\\\ \n",
        "\\vdots \\end{bmatrix}$\n",
        "\n",
        "使得  \n",
        "\n",
        "$\\pi = P\\pi$\n",
        "\n",
        "则称丌为马尔可夫裢$X = {X_{0}, X_{1},..., X(t),...}$的平稳分布。\n",
        "\n",
        "\n",
        "直观上，如果马尔可夫链的平稳分布存在，那么以该平稳分布作为初始分布，面向未来进行随机状态转移，之后任何一个时刻的状态分布都是该平稳分布。\n",
        "\n",
        "**引理19．1**\n",
        "\n",
        "给定一个马尔可夫链$X = {X_{0}, X_{1},..., X(t),...}$, 状态空间为$S$, 移概率矩阵为$P=(p_{ij})$， 则分布 $\\pi=(\\pi_{1}, \\pi_{2},...)^{T}$ 为 $X$ 的平稳分布的充要条件是$\\pi=(\\pi_{1}, \\pi_{2},...)^{T}$是下列方程组的解：\n",
        "\n",
        "$x_{i} = \\sum_{j}p_{ij}x_{j}, i=1,2,...$  \n",
        "\n",
        "$x_{i} >= 0, i = 1,2,...$  \n",
        "\n",
        "$\\sum_{i}x_{i} = 1$  \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HqSHX7PSwOjP",
        "colab_type": "text"
      },
      "source": [
        "### 吉布斯采样\n",
        "\n",
        "输入： 目标概率分布的密度函数$p(x)$, 函数$f(x)$;\n",
        "\n",
        "输出： $p(x)$的随机样本 $x_{m+1}, x_{m+2}, ..., x_{n}$，函数样本均值 $f_{mn}$;\n",
        "\n",
        "参数： 收敛步数$m$, 迭代步数 $n$.\n",
        "\n",
        "\n",
        "1. 初始化。给出初始样本 $x^{0} = $($x^{0}_{1}, x^{0}_{2},..., x^{0}_{k}$)$^{T}$.\n",
        "\n",
        "2. 对$i$循环执行  \n",
        " 设第$i-1$次迭代结束前的样本为$x^{i-1} = $($x^{i-1}_{1}, x^{i-1}_{2},..., x^{i-1}_{k}$)$^{T}$，则第$i$次迭代进行如下几步操作：  \n",
        "\n",
        " + (1)由满条件分布 $p(x_{1}|x^{i-1}_{2},...,x^{i-1}_{k})$ 抽取 $x^{i}_{1}$  \n",
        " \n",
        " + ...\n",
        " \n",
        " + (j)由满条件分布 $p(x_{j}|x^{i}_{1},...,x^{i}_{j-1}, x^{i-1}_{j+1},..., x^{i-1}_{k})$ 抽取 $x^{i}_{j}$   \n",
        " \n",
        " + (k)由满条件分布 $p(x_{k}|x^{i}_{1},...,x^{i}_{k})$ 抽取 $x^{i}_{k}$ \n",
        " \n",
        "得到第 $i$ 次迭代值 $x^{(i)} = (x^{(i)}_{1}, x^{(i)}_{2},..., x^{(i)}_{k})^{T}$.\n",
        "\n",
        "\n",
        " 3. 得到样本集合\n",
        "    \n",
        " {$x^{(m+1)}, x^{(m+2)},..., x^{(n)}$}\n",
        " \n",
        " 4. 计算\n",
        " \n",
        " $f_{mn} = \\frac{1}{n-m}\\sum_{i=m+1}^{n}f(x^{(i)})$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5ZeiXcVWBQZb",
        "colab_type": "text"
      },
      "source": [
        "--------------------------------------------------------------------------------------------------------------------------------"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PIcMXLgQBSio",
        "colab_type": "text"
      },
      "source": [
        "#### 网络资源：\n",
        "\n",
        "LDA-math-MCMC 和 Gibbs Sampling： https://cosx.org/2013/01/lda-math-mcmc-and-gibbs-sampling  \n",
        "\n",
        "MCMC蒙特卡罗方法： https://www.cnblogs.com/pinard/p/6625739.html"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kIIlKmr0I8d_",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        },
        "outputId": "265030d7-7e28-443c-eed0-e63ee8b21a73"
      },
      "source": [
        "import random\n",
        "import math\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "import numpy as np\n",
        "\n",
        "transfer_matrix = np.array([[0.6,0.2,0.2],[0.3,0.4,0.3],[0,0.3,0.7]],dtype='float32')\n",
        "start_matrix = np.array([[0.5,0.3,0.2]],dtype='float32')\n",
        "\n",
        "value1 = []\n",
        "value2 = []\n",
        "value3 = []\n",
        "for i in range(30):\n",
        "    start_matrix = np.dot(start_matrix,transfer_matrix)\n",
        "    value1.append(start_matrix[0][0])\n",
        "    value2.append(start_matrix[0][1])\n",
        "    value3.append(start_matrix[0][2])\n",
        "print(start_matrix)\n"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[[0.23076934 0.30769244 0.4615386 ]]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "A2oybaKXJGqd",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 269
        },
        "outputId": "3b2abcf8-71f1-4953-8536-667de8a1f36e"
      },
      "source": [
        "#进行可视化\n",
        "x = np.arange(30)\n",
        "plt.plot(x,value1,label='cheerful')\n",
        "plt.plot(x,value2,label='so-so')\n",
        "plt.plot(x,value3,label='sad')\n",
        "plt.legend()\n",
        "plt.show()"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VPW9//HXZyaTBJJASEKCkLDv\nBAUbxBV3tHWBLvpzaX9Ye6/6+9XbxfZWvMVbu/1q9dbahWvr7bXWul9tKWqtimhB2ypBQSFjICBI\nAAMkBEjIOvP9/TGTEAJJBhJykpn38+E8zjnfs8znMA/fc/I9y5hzDhERSQw+rwsQEZHeo9AXEUkg\nCn0RkQSi0BcRSSAKfRGRBKLQFxFJIAp9EZEEotAXEUkgCn0RkQSS5HUB7eXk5LjRo0d7XYaISL+y\nevXqPc65oV0t1+dCf/To0RQXF3tdhohIv2JmW2NZTt07IiIJRKEvIpJAFPoiIglEoS8ikkAU+iIi\nCUShLyKSQBT6IiIJpM9dpy99X9iFaQo30RBqoDHUSEOogaZQEyEXojncTLNrJhQOHZoONxNyIULh\nEM2umbALE3ZhnHORcSLjIRc6oi3swjgczrlOh8ChYZu2Fi0/C9p2mdZ57daPThyh/Ta7au9r9NOo\nfV9eWh5XTbzqhL6HQj+BHGw6yK6Du6isr6SmsYYDTQeoaayhpqmmdXig8cBh022DvTXgw01e74oc\nJ8O8LkE6MX3odIW+xKa2qZadNTupOFgRedVGhh8f/Lh1/EDjgQ7XT7Ik0pPTSQ+kk5GcQXpyOsPT\nhzPAP4BkfzIp/hSS/cmHjaf4U0jxpxDwBQj4AyT5kkiyJPzmJ8mXhN/nj0z7/PjNT8AXwGc+fObD\nb3585sPMIm1Exv3mb20zDDNrHbYsAxw23TbIWudH12kZbzs8NGizHnbY+u3nt9++SH+l0O+HDjYd\npHRvKev2rGPdnnWsr1zP1v1H3oGdnZpNXloeBRkFFOUVkZeWR97APLIHZJMRiAR7RnIG6YF0Uvwp\nCjSRBKDQ7+OaQk1s2LuB9ZXrIyFfuY5N1ZsIuzAAuQNzKcwu5MpxV1KQUUDewDzy0vLIHZBLwB/w\nuHoR6WsU+n1QU6iJv5b/lWc3PstbO99q7UPPTMlkWs40Lii4gGnZ0yjMKWTowC4fqici0kqh34ds\n3reZP278I0s3LaWqvorcgblcM/kaTh56MoXZhYxIH6EuGBHplphC38wuBX4G+IHfOOfu7mC5zwLP\nALOcc8VmNhoIAqXRRf7hnLulu0XHk4NNB3lpy0v8YeMfWLN7DUmWxLkF5/KZCZ/hrOFn4ff5vS5R\nROJIl6FvZn5gMXAxUA6sMrOlzrmSdstlAF8F3mq3iU3OuRk9VG9ccM6xbs86nt34LC9++CIHmw8y\netBobvvEbVwx7gpyBuR4XaKIxKlYjvRPA8qcc5sBzOxJYB5Q0m657wM/Bv61RyuMM8HKIN/523cI\nVgUZkDSAuaPm8pkJn2Fm7kx13YjICRdL6I8AtrWZLgdmt13AzE4FCpxzL5hZ+9AfY2bvAvuBRc65\nld0puL9yzvH4B4/zk+KfMCR1CHeefiefGvMp0pPTvS5NRBJIt0/kmpkPuA+44SizdwIjnXOVZvYJ\nYImZTXPO7W+3jZuAmwBGjhzZ3ZL6nOr6au782528vu11zss/j++f9X0yUzO9LktEElAsob8dKGgz\nnR9ta5EBFAKvR7snhgFLzexK51wx0ADgnFttZpuAicBhP4LrnHsQeBCgqKgorh4QUvxxMQtXLqSq\nvoqFpy3kusnXqRtHRDwTS+ivAiaY2RgiYX8NcF3LTOfcPqD1zKOZvQ58M3r1zlCgyjkXMrOxwARg\ncw/W32eFwiEefP9BfrX2VxRkFPDopx5lavZUr8sSkQTXZeg755rN7FbgJSKXbD7knFtvZt8Dip1z\nSztZfQ7wPTNrAsLALc65qp4ovC+rqK1g4cqFFFcUc8XYK/j26d8mLZDmdVkiIlhfe9xqUVGRKy4u\n7nrBPuqv2/7KojcX0RBqYNHpi7hy3JVelyQiCcDMVjvnirpaTnfk9pDGUCM/Xf1THg0+yuSsydw7\n515GDx7tdVkiIodR6PcA5xwLVy7kla2vcN3k67it6DZS/ClelyUicgSFfg94ZuMzvLL1FW77xG18\nsfCLXpcjItIh/UZuN22u3sw9b9/DmcPPZMG0BV6XIyLSKYV+NzSEGvjWim8xMDCQH579Q3ymf04R\n6dvUvdMN96++n9K9pSy+cLEekiYi/YIOTY/TyvKVPBp8lOunXM+c/DlelyMiEhOF/nHYU7eHRW8u\nYuKQiXz9E1/3uhwRkZipe+cYhV2YRW8soraplocueUiXZopIv6Ij/WP0aMmjvLnjTb4161uMyxzn\ndTkiIsdEoX8MgpVBfvrOT7mg4AKumniV1+WIiBwzhX6MDjYd5FsrvkVWShbfPfO7ejyyiPRL6tOP\n0T2r7mHr/q3819z/0g+giEi/pSP9GLy85WWe3fgsNxbeyOyTZne9gohIH6XQ78LHtR9z19/vojC7\nkC/P/LLX5YiIdItCvxNhF2bhyoWEwiHumXMPAV/A65JERLpFod+Jtz9+m9UVq/nmrG9SMKig6xVE\nRPo4hX4nlpQtISOQwRVjr/C6FBGRHqHQ78CBxgO8uvVVPjnmk6QmpXpdjohIj1Dod+ClLS9RH6pn\n/vj5XpciItJjFPodWFK2hHGDx1GYU+h1KSIiPUahfxSb921m7e61zBs/T3feikhcUegfxZ/K/oTf\n/FwxTidwRSS+KPTbCYVDPL/pec4ecbZ+DUtE4o5Cv52/7fgbu+p26QSuiMQlhX47S8qWkJmSybn5\n53pdiohIj1Pot7GvYR+vbXuNy8ZeRsCvRy6ISPxR6LfxwuYXaAo3qWtHROKWQr+NJWVLmJw1mclZ\nk70uRUTkhFDoR5VWlRKsCuooX0TimkI/6k+b/kSSL4lPjfmU16WIiJwwCn2gKdzEC5tf4PyC8xmS\nOsTrckREThiFPrCifAVV9VXMGzfP61JERE4ohT6RE7g5A3I4a8RZXpciInJCxRT6ZnapmZWaWZmZ\nLexkuc+amTOzojZtd0TXKzWzS3qi6J60p24PK8tXcsXYK0jyJXldjojICdVlypmZH1gMXAyUA6vM\nbKlzrqTdchnAV4G32rRNBa4BpgHDgWVmNtE5F+q5XeieFza/QMiFdNWOiCSEWI70TwPKnHObnXON\nwJPA0Tq/vw/8GKhv0zYPeNI51+Cc+xAoi26vT3DOsaRsCSfnnMzYzLFelyMicsLFEvojgG1tpsuj\nba3M7FSgwDn3wrGu66WSyhLKqsuYN14ncEUkMXT7RK6Z+YD7gG90Yxs3mVmxmRXv3r27uyXF7I9l\nfyTFn8KlYy7ttfcUEfFSLKG/HShoM50fbWuRARQCr5vZFuB0YGn0ZG5X6wLgnHvQOVfknCsaOnTo\nse3BcWoINfDihy9ywcgLGJQ8qFfeU0TEa7GE/ipggpmNMbNkIidml7bMdM7tc87lOOdGO+dGA/8A\nrnTOFUeXu8bMUsxsDDABeLvH9+I4vLbtNfY37tcJXBFJKF1eveOcazazW4GXAD/wkHNuvZl9Dyh2\nzi3tZN31ZvY0UAI0A1/uK1fuLClbwrC0YcweNtvrUkREek1MF6Y75/4M/Lld2793sOx57aZ/CPzw\nOOs7ISpqK/j7jr/zpcIv4ff5vS5HRKTXJOQduc9vfp6wC6trR0QSTkKGfnFFMeMzxzNy0EivSxER\n6VVxE/qNzWHWbKumYn99p8s55yipLGFq9tReqkxEpO+Im9CvPtjI/MVv8pd1H3e63O663VTVVyn0\nRSQhxU3oD81IIXNggNKKA50u90HVBwBMyZrSG2WJiPQpcRP6ZsbE3Aw2fNx56JdUlmAYk7Im9VJl\nIiJ9R9yEPsDEYelsqDiAc67DZYKVQUYNGkVaIK0XKxMR6RviK/TzMthf30zF/oYOlwlWBdW1IyIJ\nK+5CH2BDB/361fXV7KzdyZRshb6IJKaECv1gVRBAoS8iCSuuQj8rLZmc9JSuQ1/dOyKSoOIq9AEm\nDUuntKLmqPOClUGGpw1ncMrgXq5KRKRviLvQn5CbwcaKA4TDR17BE6wKqmtHRBJa3IX+pGEZHGwM\nsb267rD2msYatu7fqq4dEUlocRf6E/PSgSNP5pbuLQV0EldEElvchf6E1it4Du/XD1bqJK6ISNyF\n/qDUACcNTj3iSD9YFSRnQA5DB/bOb/CKiPRFMf1yVn8zMS/jqKGvo3yRvqepqYny8nLq6zt/LLpE\npKamkp+fTyAQOK714zT00/n75kpCYYffZ9Q317O5ejPnF5zvdWki0k55eTkZGRmMHj0aM/O6nD7N\nOUdlZSXl5eWMGTPmuLYRd907EDnSb2wOs7WyFoCNezcSciEd6Yv0QfX19WRnZyvwY2BmZGdnd+uv\norgM/UnDDj+Zq8cviPRtCvzYdfffKi5Df3zu4ZdtllSWMCh5EMPThntZloj0IzfccAPPPPPMCdn2\nypUrmTZtGjNmzKCurq7D5c477zyKi4t79L3jMvQHJicxMmtga+h/UPUBU7Kn6GhCRDwXCoV47LHH\nuOOOO1izZg0DBgzo1fePy9CHyMncDRUHaAo3sWHvBqZm6TdxRaRjjzzyCCeffDKnnHIKX/jCFwBY\nsWIFZ555JmPHjj3sqP/ee+9l1qxZnHzyyXznO99pbX/00Uc57bTTmDFjBjfffDOhUAiA9PR0vvGN\nb3DKKafwox/9iKeffpo777yT66+/ntdff53LL7+8dRu33norDz/88AnbzzgO/Qw2766ltLKMpnAT\nk7Mme12SiPRR69ev5wc/+AHLly9n7dq1/OxnPwNg586dvPHGGzz//PMsXLgQgJdffpmNGzfy9ttv\ns2bNGlavXs2KFSsIBoM89dRTvPnmm6xZswa/389jjz0GQG1tLbNnz2bt2rUsWrSIK6+8knvvvbd1\nfm+Ky0s2IRL6zWHHGx+tAXQSV6Q/+O5z6ynZsb9Htzl1+CC+c8W0TpdZvnw5V111FTk5OQBkZWUB\nMH/+fHw+H1OnTqWiogKIhP7LL7/MzJkzAaipqWHjxo289957rF69mlmzZgFQV1dHbm4uAH6/n89+\n9rM9ul/HK65DH+CdivUMTBrIqEGjPK5IRPqblJSU1vGW3952znHHHXdw8803H7bsL37xCxYsWMCP\nfvSjI7aTmpqK3+8/6nskJSURDodbp0/0TWpxG/pjh6bh9xmbqkuZnDUZn8VtT5ZI3OjqiPxEueCC\nC/j0pz/NbbfdRnZ2NlVVVR0ue8kll7T2x6enp7N9+3YCgQAXXngh8+bN4+tf/zq5ublUVVVx4MAB\nRo3q/IBz1KhRlJSU0NDQQF1dHa+++ipnn312T+9iq7gN/dSAn5HZqexp+pCLsvrGn1Ui0jdNmzaN\nb3/725x77rn4/f7WrpujmTt3LsFgkDPOOAOInKR99NFHmTp1Kj/4wQ+YO3cu4XCYQCDA4sWLuwz9\ngoICrr76agoLCxkzZkyn790TrOVPlr6iqKjI9dR1qQt+/xzvhP+N75/1feaPn98j2xSRnhUMBpky\nRefcjsXR/s3MbLVzrqirdeO6zyN90C4Axg2a5HElIiJ9Q1yHPsnluHAS4cZcrysREekT4jr09zZ/\nSLhhGJt3d3ybs4hIIonb0HfOsaVmA65h+BG/oiUikqhiCn0zu9TMSs2szMwWHmX+LWb2vpmtMbM3\nzGxqtH20mdVF29eY2a96egc6sr1mOwcaD5ATGMuGjw90vYKISALo8pJNM/MDi4GLgXJglZktdc6V\ntFnscefcr6LLXwncB1wanbfJOTejZ8vu2gdVHwAwdvBkSssV+iIiENuR/mlAmXNus3OuEXgSmNd2\nAedc2/um0wDPrwMtqSzBb35mDptM+d46ahuavS5JRMRzsYT+CGBbm+nyaNthzOzLZrYJuAf4SptZ\nY8zsXTP7q5md061qj0GwKsi4zHFMGZYNwMZd6tcXEemxE7nOucXOuXHA7cCiaPNOYKRzbiZwG/C4\nmQ1qv66Z3WRmxWZWvHv37p6ohZLKEiZnTWZSXsuvaKmLR0SOrra2lssuu4xTTjmFwsJCnnrqKV59\n9VVmzpzJ9OnTufHGG2loaDhivZ07dzJnzhxmzJhBYWEhK1euBOCJJ55g+vTpFBYWcvvtt/f27nQq\nltDfDhS0mc6PtnXkSWA+gHOuwTlXGR1fDWwCJrZfwTn3oHOuyDlXNHTo0Fhr79Duut1U1VcxNXsq\nBVkDSUny6WSuiHToL3/5C8OHD2ft2rWsW7eOSy+9lBtuuIGnnnqK999/n+bmZh544IEj1nv88ce5\n5JJLWLNmDWvXrmXGjBns2LGD22+/neXLl7NmzRpWrVrFkiVLPNiro4vl2TurgAlmNoZI2F8DXNd2\nATOb4JzbGJ28DNgYbR8KVDnnQmY2FpgAbO6p4jvSchJ3StYU/D5jQl46G9S9I9L3vbgQPn6/Z7c5\nbDp88u5OF5k+fTrf+MY3uP3227n88ssZNGgQY8aMYeLEyDHqggULWLx4MV/72tcOW2/WrFnceOON\nNDU1MX/+fGbMmMHy5cs577zzaDmAvf7661mxYgXz5/eNR8F0eaTvnGsGbgVeAoLA08659Wb2veiV\nOgC3mtl6M1tDpBtnQbR9DvBetP0Z4BbnXMePr+shJZUlGMakrMjjFybmZehIX0Q6NHHiRN555x2m\nT5/OokWLOjwyf+utt5gxYwYzZsxg6dKlzJkzhxUrVjBixAhuuOEGHnnkkV6u/NjF9JRN59yfgT+3\na/v3NuNf7WC9Z4Fnu1Pg8QhWBhk1aBRpgTQgEvp/eGc7++qaGDwg0NvliEisujgiP1F27NhBVlYW\nn//858nMzOSXv/wlW7ZsoaysjPHjx/P73/+ec889l9mzZ7NmzZrW9bZu3Up+fj7//M//TENDA++8\n8w633347X/nKV9izZw9DhgzhiSee4F/+5V882a+jictHKwergswYeujWgJaTuRsrDlA0OsurskSk\nj3r//ff513/9V3w+H4FAgAceeIB9+/Zx1VVX0dzczKxZs7jllluOWO/111/n3nvvJRAIkJ6eziOP\nPMJJJ53E3Xffzfnnn49zjssuu4x58+Yd5V29EXehX11fzc7anVw7+drWtgl56QCUKvRF5CguueQS\nLrnkkiPa33333U7XW7BgAQsWLDii/dprr+Xaa689yhrei7tn7wSrgsDhv4k7InMAacl+NuoZPCKS\n4OI39LMOhb6ZMSEvg1KdzBWRBBd/oV8ZZHjacAanDD6sfVJeBht3KfRFJLHFX+hXBQ/r2mkxIS+d\nPTWNVNYceVediEiiiKvQr2msYev+rYd17bSYNKzlcQzq1xeRxBVXoV+6txTgqEf6egaPiEichX6w\n8siTuC2GZqQweECAUoW+iPSQLVu2UFhY6HUZxyS+Qr8qSM6AHIYOPPKhbWYWOZmr0BeRBBZXN2cF\nq4JHPcpvMSEvnefW7sA5h5n1YmUi0pfV1tZy9dVXU15eTigU4s4776S0tJTnnnuOuro6zjzzTH79\n619jZqxevZobb7wRgLlz53pc+bGLmyP9+uZ6NldvPmp/fotJwzLYX9/MrgO6gkdEDjnao5VvvfVW\nVq1axbp166irq+P5558H4Itf/CK/+MUvWLt2rcdVH5+4OdKvaaph7ui5FOUVdbjMhNzIydzSjw+Q\nNyi1t0oTkRj9+O0ftz4avadMzprM7ad1/kMm7R+tfM455/Dss89yzz33cPDgQaqqqpg2bRrnnHMO\n1dXVzJkzB4AvfOELvPjiiz1a74kWN6GfMyCHe+bc0+kyE6PP4NlQcYA5E7v/Yy0iEh9aHq385z//\nmUWLFnHhhReyePFiiouLKSgo4K677qK+vt7rMntE3IR+LLLTU8hJT9FlmyJ9VFdH5CdK+0cr/+Y3\nvwEgJyeHmpoannnmGT73uc+RmZlJZmYmb7zxBmeffTaPPfaYJ/V2R0KFPkSO9kt1g5aItHG0Rysv\nWbKEwsJChg0bxqxZs1qX/e1vf8uNN96ImfXLE7nmnPO6hsMUFRW54uLiE7b9u5au53+Kt/H+XZfg\n8+kKHhGvBYNBpkzp+AIMOdLR/s3MbLVzruOTmlFxc/VOrCbmZVDbGGJ7dZ3XpYiI9LqEC/1JwyIn\nc/XETRFJRAkX+uNbL9tUv76IJJ6EC/3BAwKcNDhVj2MQ6UP62rnFvqy7/1YJF/oQ6dfXg9dE+obU\n1FQqKysV/DFwzlFZWUlq6vHfXJpwl2wCFI4YxK//upn99U0MSg14XY5IQsvPz6e8vJzdu3d7XUq/\nkJqaSn5+/nGvn5Chf/6kXBa/tom/lu7milOGe12OSEILBAKMGTPG6zISRkJ278wcOYSstGSWBSu8\nLkVEpFclZOj7fcb5k3J5vXQ3zaGw1+WIiPSahAx9gIun5rKvronirXu9LkVEpNckbOifM2EoyX4f\ny0rUxSMiiSNhQz8tJYkzxmWzLFihS8VEJGEkbOgDXDQlly2VB9m0u9brUkREekVCXrLZ4sIpedz5\np/W8GqxgfG661+V0X6gZmuuguQGaosOjTYcaI8uGmyHcFBl2NO1C4MIQDoFzkfGWttb28KF5uBjG\nOTTdonXcHT5+xHximNdu/mHNsf5Vp7/+xAM5E+GTPz6hb5HQoT88cwBTTxrEsmAFN587zutyDmmq\nh/3bYd82qNkF9fugvjo67OQVbu7ZOswH5gef/9C4+cCsTVtLu0XGsei4Rcd9Rx+HTsZp195aULv6\nOpl3xPzDZsS4/3r0tvSyxhPf65DQoQ+RLp5fvlbG3tpGhqQl986b1u6BvVsjob6v/FDA7yuPvGo7\nuDMxMBBSBx96pQ2F7PGR8ZRBkfmBVEhKhaQUSBoQGQaiw6ToPH8y+APgS4q82o63TJsffAnd+ycS\nlxT6U/P4+fIyXivdxWdOPf5bmzvU3AA734PyVZHX9mKo/ujwZQJpkFkAg/PhpFMiw8HR6fRhMCAz\nEupJvfSlJCJxK6bQN7NLgZ8BfuA3zrm7282/BfgyEAJqgJuccyXReXcAX4rO+4pz7qWeK7/7CocP\nJjcjhWXBiu6HvnNQvRXKi6OvVfDxe5E+dIgE+YhPwGk3RY7QB+dHXqmZ6koQkV7RZeibmR9YDFwM\nlAOrzGxpS6hHPe6c+1V0+SuB+4BLzWwqcA0wDRgOLDOzic65UA/vx3Hz+YwLp+Tx3NodNDSHSEny\nH9sGnIsE/DsPw4aXDnXNJA2AEafC6f8H8mfBiCIYdFKP1y8icixiOdI/DShzzm0GMLMngXlAa+g7\n5/a3WT6NQ5c+zAOedM41AB+aWVl0e3/vgdp7zEVTcnni7Y94a3MVcyYOjW2lur3w3tOw+newa32k\ni2byZTBydiTkc6eBP+F7z0Skj4kllUYA29pMlwOz2y9kZl8GbgOSgQvarPuPduuOOK5KT6CzxueQ\nGvDxarCi89B3Dj76B6x+GEqWQHM9DJ8Jl98P0z8HKRm9VrOIyPHosUNR59xiYLGZXQcsAhbEuq6Z\n3QTcBDBy5MieKilmqQE/Z48fyrLgLu660mHt+9cPVsHaJyJH9XtKITkDZlwPn1gQOfEqItJPxBL6\n24GCNtP50baOPAk8cCzrOuceBB4EKCoq8uSumIun5rIsWMEHHx9gykmDIo37d8Cy78L6P0ROxubP\ngnmLYdqnITnNizJFRLolltBfBUwwszFEAvsa4Lq2C5jZBOfcxujkZUDL+FLgcTO7j8iJ3AnA2z1R\neE87f3IuAMtKKpgyLAPWPAZ/+bdI2H/ii5Gj+rxpHlcpItI9XYa+c67ZzG4FXiJyyeZDzrn1ZvY9\noNg5txS41cwuApqAvUS7dqLLPU3kpG8z8OW+dOVOW7kZqcwoyGTN+vWwcyGULYORZ8K8X0J2H7pb\nV0SkG6yvPWGyqKjIFRcX9/4bO8fyx+9l1ob7SA+AXfw9mPVPuitVRPoFM1vtnCvqajldUwiRRyI8\n9xUu2Pw6fwtPpfLs/+CK2Wd5XZWISI9L7NAPh2H1Q/DKdwBwl93Ht5aNYvLWFK7wuDQRkRMhcUO/\n6kNY+i+wZSWMPR+u/DmWOZILd6zjqeJt1DeFSA0c4925IiJ9XOJ1WDsHq/4bHjgTdq6FK34OX/gj\nZEbuD7hoah71TWHeLNvjcaEiIj0vsUK/bi88/b/hhdtg5Onwf/8euRSzzc1Ys8dkk56SxLKgfjtX\nROJP4nTvbHsbnvkSHNgBF38fzrj1qFfmJCf5OHfiUF4N7iIcdvh8evqliMSP+D/SD4dh5U/goUsj\nR/Q3vgxnfaXTSzEvnJLLrgMNrNuxrxcLFRE58eL7SP9ABfzxJtj8Okz7DFxxf+RXprpw/qRcfBa5\nO/fk/MwTX6eISC+J3yP9slfhV2fBR29FTtZ+7qGYAh9gSFoyRaOyWBbcdYKLFBHpXfEX+qGmyHX3\nj34GBubATa8dcbI2FhdNzaVk5362V9edoEJFRHpffIX+3i2Rvvs37488JO2m1yB3ynFt6sIpeQAs\n11U8IhJH4if095TBr+bAno1w1cOR/vvAgOPe3Lih6YzJSeMVdfGISByJn9DPGgtFX4RbVkSed98D\nLpqSyz82VVLT0Nwj2xMR8Vr8hL7PBxd/F4aM7rFNXjglj8ZQmFfVxSMicSJ+Qv8EKBo1hLFD0/jZ\nqxtpCoW9LkdEpNsU+p1I8vv4t09OYfPuWp54+yOvyxER6TaFfhcunJLLGWOz+ekrG9hX1+R1OSIi\n3aLQ74KZ8e3LplBd18R/vlbmdTkiIt2i0I9B4YjBfPbUfH775ha2VR30uhwRkeOm0I/RN+dOwu8z\n7v7LB16XIiJy3BT6MRo2OJWb5ozlhfd2snprldfliIgcF4X+Mbj53LHkZqTwgxeCOOe8LkdE5Jgp\n9I/BwOQkvnnJJN79qJrn39vpdTkiIsdMoX+MPntqPlNOGsTdL35AfVPI63JERI6JQv8Y+X3Gosum\nsL26jof/tsXrckREjolC/zicNT6HCyfnsnh5GZU1DV6XIyISM4X+cbrjU1M42BTi/mUbvS5FRCRm\nCv3jND43netnj+Txtz+ibNf3Gr1SAAAJ9klEQVQBr8sREYmJQr8bvnrhBAYG/Py/P+uGLRHpHxT6\n3ZCdnsKtF4xn+Qe7eGPjHq/LERHpkkK/mxacOZr8IQP4wQslhMK6YUtE+jaFfjelBvws/ORkPvj4\nAM+s3uZ1OSIinVLo94DLpp/EqSMz+dGLH7Bu+z6vyxER6ZBCvweYGff/r5mkJSdx7YP/0APZRKTP\niin0zexSMys1szIzW3iU+beZWYmZvWdmr5rZqDbzQma2Jvpa2pPF9yUjswfyP7ecQU5GCp//zdu8\nWaYTuyLS93QZ+mbmBxYDnwSmAtea2dR2i70LFDnnTgaeAe5pM6/OOTcj+rqyh+ruk4ZnDuCpm09n\nZNZAvvjwKpaVVHhdkojIYWI50j8NKHPObXbONQJPAvPaLuCce8051/KTUv8A8nu2zP4jNyOVJ286\nncnDMrjl0dU8t3aH1yWJiLSKJfRHAG0vSymPtnXkS8CLbaZTzazYzP5hZvOPtoKZ3RRdpnj37t0x\nlNS3DUlL5rF/ms2pI4fw1Sff5eliXdUjIn1Dj57INbPPA0XAvW2aRznnioDrgPvNbFz79ZxzDzrn\nipxzRUOHDu3JkjyTkRrgdzeexlnjc/jWM+/x8Jsfel2SiEhMob8dKGgznR9tO4yZXQR8G7jSOdf6\n6Enn3PbocDPwOjCzG/X2KwOS/fxmQRFzp+Zx13MlLH6tzOuSRCTBxRL6q4AJZjbGzJKBa4DDrsIx\ns5nAr4kE/q427UPMLCU6ngOcBZT0VPH9QUqSn8XXn8q8GcO596VS7n3pA/3Uooh4JqmrBZxzzWZ2\nK/AS4Acecs6tN7PvAcXOuaVEunPSgf8xM4CPolfqTAF+bWZhIl8wdzvnEir0AQJ+H/ddPYOByX4W\nv7aJ2oYQ/375VHw+87o0EUkw1teOOouKilxxcbHXZZwQzjm+/3yQh978kOkjBvO1iyZwweRcol+U\nIiLHzcxWR8+fdkp35PYiM+POy6fwk6tOobqukS/9rpj5i9/ktdJd6vIRkV6hI32PNIXC/OGdcn7+\nahnbq+uYOTKTr180kXMm5OjIX0SOWaxH+gp9jzU2h3lmdTm/XL6RHfvqKRo1hK9fPJEzx2Ur/EUk\nZgr9fqahOcTTxeUsXl7Gx/vrOW1MFl+/aCJnjMv2ujQR6QcU+v1UfVOIJ9/+iP98fRO7DjRwSkEm\n508ayjkTcjg5P5OAX6dhRORICv1+rr4pxONvfcSSNdt5f/s+nIP0lCROH5vNORNyOHtCDmNz0tQF\nJCKAQj+uVB9s5G+bKlm5cQ9vlO1mW1UdAMMHp3LW+MgXwFnjc8hJT/G4UhHxikI/jn1UeZCVZbt5\nY+Me/rapkn11TQAMzUhhTHYao7IHMjonjdFtxtNTurwPT0T6MYV+ggiFHeu27+PvmyvZtKuGrZUH\n+bCylt0HGg5bLic9hTE5AxmVnUb+kAEMGZhM5sAAWWnJh40PCPjVZSTSD8Ua+jr86+f8PuOUgkxO\nKcg8rL22oZmtlQfZUlkbee2pZUvlQVZs2M2udl8IbSUn+RgyMMCQgZEvg7QUP6kBPwMCfgYkR19t\np6PDlCQ/Ab8R8PtI8hmBJB8Bn49AkpHk85Hs95Hkt8jL58Nn4PMZfjP8PsPXOkRfOiInkEI/TqWl\nJDF1+CCmDh90xLzmUJjquiaqDzay92ATVbWNreN7DzaytzYyXn2wkZ37mqhrClHXGGodNjSHT2jt\nZuA3w+czDPCZYRYdRue3nweRobXZhrVpa/tF0jLaOoyudWi6Zf6RXz4dfh3F+D11LF9n/eXLr39U\n2T9MPmkQv7j2xD6IWKGfgJL8PnLSU477xG8o7KhvOvQlUN8Uor4pTFM4TFNzmOawozEUpjnkaA6F\nW8ebQmGawo5w2BEKO8IuMgy5ljYIu8PbifxHOOwiQ+dwLvIco7ADR3ToILIk0fmRedFN0NKL6Wgd\naTtofQzGoekj97ujjtBYu0iPqSO1b/W6dsj1l0L7iYIhA074eyj05Zj5fUZaShJpOjks0u/oTh8R\nkQSi0BcRSSAKfRGRBKLQFxFJIAp9EZEEotAXEUkgCn0RkQSi0BcRSSB97oFrZrYb2NqNTeQAe3qo\nnL4g3vYH4m+f4m1/IP72Kd72B47cp1HOuaFdrdTnQr+7zKw4lifN9Rfxtj8Qf/sUb/sD8bdP8bY/\ncPz7pO4dEZEEotAXEUkg8Rj6D3pdQA+Lt/2B+NuneNsfiL99irf9gePcp7jr0xcRkY7F45G+iIh0\nIG5C38wuNbNSMyszs4Ve19MTzGyLmb1vZmvMrN/9cLCZPWRmu8xsXZu2LDN7xcw2RodDvKzxWHWw\nT3eZ2fbo57TGzD7lZY3HwswKzOw1Mysxs/Vm9tVoe7/8nDrZn/78GaWa2dtmtja6T9+Nto8xs7ei\nmfeUmSXHtL146N4xMz+wAbgYKAdWAdc650o8LaybzGwLUOSc65fXF5vZHKAGeMQ5Vxhtuweocs7d\nHf1yHuKcu93LOo9FB/t0F1DjnPsPL2s7HmZ2EnCSc+4dM8sAVgPzgRvoh59TJ/tzNf33MzIgzTlX\nY2YB4A3gq8BtwB+cc0+a2a+Atc65B7raXrwc6Z8GlDnnNjvnGoEngXke15TwnHMrgKp2zfOA30XH\nf0fkf8h+o4N96recczudc+9Exw8AQWAE/fRz6mR/+i0XUROdDERfDrgAeCbaHvNnFC+hPwLY1ma6\nnH7+QUc54GUzW21mN3ldTA/Jc87tjI5/DOR5WUwPutXM3ot2//SLrpD2zGw0MBN4izj4nNrtD/Tj\nz8jM/Ga2BtgFvAJsAqqdc83RRWLOvHgJ/Xh1tnPuVOCTwJejXQtxw0X6Fvt//yI8AIwDZgA7gZ94\nW86xM7N04Fnga865/W3n9cfP6Sj7068/I+dcyDk3A8gn0rMx+Xi3FS+hvx0oaDOdH23r15xz26PD\nXcAfiXzY/V1FtN+1pf91l8f1dJtzriL6P2UY+C/62ecU7Sd+FnjMOfeHaHO//ZyOtj/9/TNq4Zyr\nBl4DzgAyzSwpOivmzIuX0F8FTIiezU4GrgGWelxTt5hZWvREFGaWBswF1nW+Vr+wFFgQHV8A/MnD\nWnpESzhGfZp+9DlFTxL+NxB0zt3XZla//Jw62p9+/hkNNbPM6PgAIhesBImE/+eii8X8GcXF1TsA\n0Uuw7gf8wEPOuR96XFK3mNlYIkf3AEnA4/1tn8zsCeA8Ik8DrAC+AywBngZGEnma6tXOuX5zYrSD\nfTqPSLeBA7YAN7fpD+/TzOxsYCXwPhCONv8bkX7wfvc5dbI/19J/P6OTiZyo9RM5UH/aOfe9aEY8\nCWQB7wKfd841dLm9eAl9ERHpWrx074iISAwU+iIiCUShLyKSQBT6IiIJRKEvIpJAFPoiIglEoS8i\nkkAU+iIiCeT/A4c3KvXb4eVeAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "z2OC3tyoJNeN",
        "colab_type": "text"
      },
      "source": [
        "可以发现，从10轮左右开始，我们的状态概率分布就不变了，一直保持在  \n",
        "[0.23076934,0.30769244,0.4615386]\n",
        "\n",
        "### https://zhuanlan.zhihu.com/p/37121528"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3Taf1Z8nJU12",
        "colab_type": "text"
      },
      "source": [
        "#### M-H采样python实现  \n",
        "https://zhuanlan.zhihu.com/p/37121528\n",
        "\n",
        "假设目标平稳分布是一个均值3，标准差2的正态分布，而选择的马尔可夫链状态转移矩阵 $Q(i,j)$ 的条件转移概率是以 $i$ 为均值,方差1的正态分布在位置 $j$ 的值。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "eUoFG0tYJx9e",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 323
        },
        "outputId": "b404de1a-535b-4772-f90d-91f1f17ebe1e"
      },
      "source": [
        "from scipy.stats import norm\n",
        "\n",
        "def norm_dist_prob(theta):\n",
        "    y = norm.pdf(theta, loc=3, scale=2)\n",
        "    return y\n",
        "\n",
        "T = 5000\n",
        "pi = [0 for i in range(T)]\n",
        "sigma = 1\n",
        "t = 0\n",
        "while t < T-1:\n",
        "    t = t + 1\n",
        "    pi_star = norm.rvs(loc=pi[t - 1], scale=sigma, size=1, random_state=None)   #状态转移进行随机抽样\n",
        "    alpha = min(1, (norm_dist_prob(pi_star[0]) / norm_dist_prob(pi[t - 1])))   #alpha值\n",
        "\n",
        "    u = random.uniform(0, 1)\n",
        "    if u < alpha:\n",
        "        pi[t] = pi_star[0]\n",
        "    else:\n",
        "        pi[t] = pi[t - 1]\n",
        "\n",
        "\n",
        "plt.scatter(pi, norm.pdf(pi, loc=3, scale=2),label='Target Distribution')\n",
        "num_bins = 50\n",
        "plt.hist(pi, num_bins, normed=1, facecolor='red', alpha=0.7,label='Samples Distribution')\n",
        "plt.legend()\n",
        "plt.show()"
      ],
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.6/dist-packages/matplotlib/axes/_axes.py:6521: MatplotlibDeprecationWarning: \n",
            "The 'normed' kwarg was deprecated in Matplotlib 2.1 and will be removed in 3.1. Use 'density' instead.\n",
            "  alternative=\"'density'\", removal=\"3.1\")\n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl41OW99/H3l7AEkR2EshlUSFib\nYAQsyqKgVEWo1YqVR63bUY96aJU+sXi5UGk5pRfaHqEt5VisRaVSzQOVloKgrQvKFmUXUISgAiIg\nyBq4nz9mkoYwydxJZp/P67pykfn97t9vvgmT79xzr+acQ0RE0kOdeAcgIiKxo6QvIpJGlPRFRNKI\nkr6ISBpR0hcRSSNK+iIiaURJX0QkjSjpi4ikESV9EZE0UjfeAVTUqlUrl5WVFe8wRESSyooVK75w\nzrUOVy7hkn5WVhbLly+PdxgiIknFzD7xKafmHRGRNKKkLyKSRpT0RUTSSMK16Ydy/PhxiouLOXLk\nSLxDkSSQmZlJhw4dqFevXrxDEUk4SZH0i4uLady4MVlZWZhZvMORBOacY8+ePRQXF9O5c+d4hyOS\ncJKieefIkSO0bNlSCV/CMjNatmypT4UilUiKpA8o4Ys3vVZEKpc0SV9ERGovKdr0TzNiRGTvN29e\n2CITJ07k+eefJyMjgzp16vC73/2Ofv36RTaOcgYPHswvf/lL8vPza3yPmTNnMm7cODp06MDBgwc5\n55xzePTRR/nWt74FwCOPPMLAgQMZOnRoyOsLCwvp2rUr3bt3D3n+t7/9LWeccQY33XRTtePdt28f\nzz//PPfccw8An376Kffffz9z5sypwU8qIr6SM+nH2DvvvMNf//pXVq5cSYMGDfjiiy84duxYvMPy\ncv311/P0008DsGTJEq655hqWLFlCt27dmDBhQpXXFhYWctVVV4VM+iUlJdx11101jmvfvn1Mmzat\nLOm3a9dOCV8kBtS84+Gzzz6jVatWNGjQAIBWrVrRrl07ACZMmMAFF1xAz549ufPOO3HOAYGa+g9/\n+EPy8/Pp1q0by5Yt45prrqFLly48/PDDAGzdupWcnBxuvPFGunXrxrXXXsuhQ4dOe/5//OMfXHjh\nhfTp04frrruOgwcPAlBQUED37t3p3bs3Dz74YNifY8iQIdx5551Mnz4dgFtuuaUs0Va819tvv83c\nuXMZN24cubm5bNmyhcGDBzN27Fjy8/P51a9+xWOPPcYvf/nLsvs/99xz5Obm0rNnT9577z2A08r0\n7NmTrVu3UlBQwJYtW8jNzWXcuHFs3bqVnj17AoGO+x/84Af06tWLvLw8lixZAgQ+uVxzzTUMHz6c\nLl268OMf/9j3v1BEgpT0PVx22WVs376drl27cs899/DGG2+Unbv33ntZtmwZa9as4fDhw/z1r38t\nO1e/fn2WL1/OXXfdxciRI5k6dSpr1qxh5syZ7NmzB4CNGzdyzz33sH79epo0acK0adNOee4vvviC\nJ554gkWLFrFy5Ury8/OZMmUKe/bs4ZVXXmHt2rV88MEHZW8k4fTp04cNGzaccizUvb71rW9x9dVX\nM3nyZIqKijj33HMBOHbsGMuXL+eBBx447d6HDh2iqKiIadOmceutt1YZx6RJkzj33HMpKipi8uTJ\np5ybOnUqZsbq1at54YUXuPnmm8tG4xQVFTF79mxWr17N7Nmz2b59u9fPLSIBSvoezjzzTFasWMH0\n6dNp3bo1119/PTNnzgQCTSb9+vWjV69eLF68mLVr15Zdd/XVVwPQq1cvevTowTe+8Q0aNGjAOeec\nU5asOnbsyIABAwAYM2YMb7755inPvXTpUtatW8eAAQPIzc3l2Wef5ZNPPqFp06ZkZmZy22238fLL\nL3PGGWd4/Syln0TKq869rr/++krP3XDDDQAMHDiQr776in379nnFVNGbb77JmDFjAMjJyeHss8/m\nww8/BODSSy8ti7d79+588onXGlMiEqQ2fU8ZGRkMHjyYwYMH06tXL5599llGjx7NPffcw/Lly+nY\nsSOPPfbYKePDS5uD6tSpU/Z96eOSkhLg9OGFFR875xg2bBgvvPDCaTG99957vPbaa8yZM4enn36a\nxYsXh/05Vq1aRbdu3U45VrduXe97NWrUqNJ7h/pZ6taty8mTJ8uO1Xb8fPnfY0ZGRtnvMWFVNujA\nY/CASDSopu9h48aNbNq0qexxUVERZ599dlkCa9WqFQcPHqxRR+S2bdt45513AHj++ee56KKLTjnf\nv39/3nrrLTZv3gzA119/zYcffsjBgwfZv38/V1xxBU8++STvv/9+2Od64403mD59Onfccccpxyu7\nV+PGjTlw4ID3zzJ79mwgUFNv2rQpTZs2JSsri5UrVwKwcuVKPv7447D3vvjii5k1axYAH374Idu2\nbSM7O9s7DhGpXHLW9GNcSzp48CD33Xcf+/bto27dupx33nlMnz6dZs2acccdd9CzZ0/atm3LBRdc\nUO17Z2dnM3XqVG699Va6d+/O3Xfffcr51q1bM3PmTG644QaOHj0KwBNPPEHjxo0ZOXIkR44cwTnH\nlClTQt5/9uzZvPnmmxw6dIjOnTvzl7/85bSa/oEDB0Lea/To0dxxxx38+te/9npDy8zMJC8vj+PH\nj/PMM88A8N3vfpc//vGP9OjRg379+tG1a1cAWrZsyYABA+jZsyff/va3+c///M+y+9xzzz3cfffd\n9OrVi7p16zJz5sxTavgiUnMWqo33tEJmw4FfARnADOfcpArnfwTcDpQAu4FbnXOfBM/dDJT2Mj7h\nnHu2qufKz893FTdRWb9+/WmJKhVs3bqVq666ijVr1sQ7lJSTMK8ZNe9IjJjZCudc2IkyYWv6ZpYB\nTAWGAcXAMjOb65xbV67YKiDfOXfIzO4GfgFcb2YtgEeBfMABK4LX7q3+jySSmLIKXq303Iz1OwEY\n2q1NrMIRqZJP805fYLNz7iMAM3sRGAmUJX3n3JJy5ZcCY4LfXw4sdM59Gbx2ITAcOL1XMg1lZWWp\nlp+kqkr0oSwKJn/QG4DEl09Hbnug/GDo4uCxytwG/K2G14oktKyCV6ud8CtatH5nre8hUlMR7cg1\nszEEmnIGVfO6O4E7ATp16hTJkEQiImf8fI6cCN//VR1ZBa8ypn8nnhjVK6L3FamKT01/B9Cx3OMO\nwWOnMLOhwHjgaufc0epc65yb7pzLd87lt27d2jd2kZjIKng14gm/1J+WblOtX2LKJ+kvA7qYWWcz\nqw+MBuaWL2BmecDvCCT8XeVOLQAuM7PmZtYcuCx4TCQpxCohK/FLrIRN+s65EuBeAsl6PfBn59xa\nM5tgZlcHi00GzgReMrMiM5sbvPZL4KcE3jiWARNKO3WTyZ49e8jNzSU3N5e2bdvSvn37ssfRWm1z\n5cqV/P3vfw95btGiRTRt2pS8vDy6du3KoEGDmD9/ftn5qVOnlk1uCmXx4sUsXbq00vOvvPJK2Xo4\nY8aMobCw0DvukydPMmnSv0f0njhxgosvvtj7+kQSzUQ8Y87jp30tOq9v1J5PpJRXm75zbj4wv8Kx\nR8p9H3pB9sC5Z4BnahpgImjZsiVFRUVAYNXIM88802tVy1InTpwgIyOjWs+5cuVK1qxZw/Dhw0Oe\nHzJkSFkyXrlyJd/5znf44x//yKBBg06Z6BTK4sWLadWqFf379z/tXElJCd/5zneqFWt5pUm/oKAA\nCCyV8K9//avG94uXeNW8swpeZeukK+Py3JIeUnIZhsJVOxgwaTGdC15lwKTFFK46rRshYkaMGMH5\n559Pjx49mDFjBhBInM2aNWPs2LH07t2b9957j7lz55Kdnc3555/Pfffdx6hRo4DAbN9bbrmFvn37\nkpeXx7x58zh8+DATJkxg1qxZ5Obmhp0N26dPH8aPH1+2bv7DDz/MU089BcCTTz5ZtmTymDFj2LJl\nCzNmzGDy5Mnk5uby9ttvM2bMGO6++2769u3LT37yE2bMmMHYsWPL7r9gwQLOP/98unbtyt/+FhiY\nVbHM8OHDefPNNykoKODAgQPk5uZy0003lf0uIPCG8KMf/YiePXvSq1evsp9r0aJFXHrppVxzzTVk\nZ2dz0003ReK/psZqkvAzM4ytk6487StWzy/iKzmXYahC4aodPPTyag4fPwHAjn2Heejl1QCMyov8\naNFnn32WFi1acOjQIfLz8/nud79L48aN2b9/PwMHDuSpp57i0KFDdO3albfeeotOnTrxve99r+z6\nCRMmMHz4cGbOnMnevXvp168fH3zwAY888ghr1qwpS97h9OnTh//5n/857fgvfvELPvnkE+rXr8++\nffto1qwZt99+O61atSpL2tOmTeOzzz5j6dKl1KlTp+zNq9T27dtZtmwZmzZtYujQoWXrAIUyadIk\nZsyYUfbJqPyCaC+99BLr16/n/fffZ/fu3VxwwQUMHDgQCHxaWbt2LW3atKF///4sXbo05CeRaKtJ\nwq0quZeOyS8/Tt83DtX4JRpSrqY/ecHGsoRf6vDxE0xesDEqz/fkk0/yzW9+kwsvvJDi4mK2bNkC\nBNbSL20mWbduHdnZ2Zx99tmYWdkSxBDYIGXixInk5uYyZMgQjhw5wrZt26odR2XLafTo0YMxY8Yw\na9Ys6tWrV+n11113HXXqhH45fO9736NOnTpkZ2fTsWPHUxafq44333yTG264gYyMDNq2bctFF11E\n6ZIb/fv3p127dmRkZJCbm8vWrVtr9By10W/iwmqVr05tfmi3NpqUJQkh5Wr6n+47XK3jtbFo0SL+\n+c9/snTpUho2bMhFF11UtvJmw4YNT1tqOBTnHIWFhWWblJT65z//Wa1YQi2ZDIGmmTfeeIO5c+fy\ns5/9jA8++CDk9VoyGXYe8O+UPy3Ze+7bPLRbG+9av2r7Eg0pV9Nv16xhtY7Xxv79+2nRogUNGzZk\n7dq1LFu2LGS57t27s3HjRrZv345zrmwJYoDLL7/8lGaZVatWAdVb1rioqIif/exnp3XgnjhxguLi\nYi655BJ+8Ytf8MUXX3Do0KFqL5n80ksv4Zzjww8/ZPv27XTp0oWsrCxWrVqFc46tW7eyYsUKILA2\nPxAyaV988cW8+OKLnDx5kp07d/LWW2/VauP3SKpOs05tE3F1avxq35dIS7mkP+7ybBrWO3WkTMN6\nGYy7PPLrsV955ZUcOnSI7t278/DDD9OvX7+Q5c444wyefvpphg4dSn5+Ps2aNaNp06YAPProo3z9\n9ddlu2s99thjAFxyySW8//775OXlhezIXbJkCXl5eWRnZ3P//fczbdo0Bg06dSJ0SUkJ3//+9+nd\nuzd9+vThwQcfLFuS+c9//jN5eXm8/fbbYX/O9u3bk5+fz4gRI5g+fTr169dn0KBBtG/fnm7duvHA\nAw+Qm5tbVv62226jd+/ep3XIXnvtteTk5NC7d2+GDh3KlClTOOuss8I+f7TFMuGXUuKXePFaWjmW\nIrG0cuGqHUxesJFP9x2mXbOGjLs8OyqduNVx8OBBzjzzTJxz/Md//Ae9evXivvvui2tMqcz3NTNs\nyuts2vW11z2rTPiezTsV+Tb1DN38Xo3uL+kjYksrJ6NRee3jnuQr+s1vfsOsWbM4evQo+fn5p+1e\nJfFR7YRfw+ReGd82frXvS6SkXPNOoho3bhxFRUWsX7+e5557jszMzHiHlPZ8m00GnNsiqnH0bNfU\nq5yaeSQSkqam75zzGg0j4tNkGW545ow5j5d9P7RbmwqrTUVW26aZrPl0f/SeQKScpKjpZ2ZmsmfP\nHq8/Zklvzjn27NkT9pOU7/DMWI2t930e1faltpKipt+hQweKi4vZvXt3vEORJJCZmUmHDh0qPe+b\nOGM9mcq3fb9zwat8rPZ9qaGkSPr16tWjc+fO8Q5D0kijBon7p6HPu1IbSdG8IxIpvrX8C89pGeVI\nQvP9dNFZzTxSQ0r6kjZ8V1uN99BIn8Sv2r7UlJK+pI2xs4viHUJEqVNXakJJX9JC70dD70JWUbxr\n+aV8m3mGTXk9uoFIyknc3iqRCPrq6ImQx8uPx2/RqD6s/W2sQgqrQ/OGFO+tenVY3xnFIqWU9CXl\nZRW8ekpyr0yfTs1jEI2/nLZNwiZ9gBt//w6z7rgwBhFJKlDzjgixH5Pvyyeut7Z8GYNIJFUo6UtK\nS5fOTt8+CxElfUl7iVrLL+UTX2V9FiIVKelLykqXWn6pdPt5pWaU9CWtJXotv1SiDCWV5KfRO5Ja\ngpucvLFxFzNOVj1vNRUX6tZmKxKOavqSko6HSfgAlyZJLb+UkrlEgpK+pByf5Ynr1UnOer5P1Df+\n/p2oxyHJS0lf0tKg7LPiHUKN+Kyjr3H7UhUlfUkp/9oUfqOdRF4rP1LCbQcp6Sv1X/2SeoKdtaEc\nLTkZ9vJ4rZUfKVsnXRl2eKbvdpCSflTTl5Txzkd7wpbJsORsy68JrcApoSjpS8r4+mhJ2DJDcpKz\nLb8in5E8WoFTQlHzjiSuKppxKnpj466wZZJ1xE5tFK7awai89vEOQxKIavqSEnzG5SfriJ3K+NT2\nU223MKk9JX2RJJaZkX6fXqR2vJK+mQ03s41mttnMCkKcH2hmK82sxMyurXDuhJkVBb/mRipwkVI+\nk7GSZY2d6tow8YqwZc57SAuxyb+FbdM3swxgKjAMKAaWmdlc59y6csW2AbcAD4a4xWHnXG4EYhWR\nGigJ3/IlacSnpt8X2Oyc+8g5dwx4ERhZvoBzbqtz7gMg/CBpkQhK51p+KZ+2/c5adlmCfJJ+e2B7\nucfFwWO+Ms1suZktNbNRoQqY2Z3BMst37w4/o1JEqkeVfSkVi47cs51z+cD3gafM7NyKBZxz051z\n+c65/NatW8cgJEkFquX/m09tXwuxCfiN098BdCz3uEPwmBfn3I7gvx+Z2etAHrClGjGKSGVzFubN\n876FFmIT8KvpLwO6mFlnM6sPjAa8RuGYWXMzaxD8vhUwAFhX9VUi4flMxurZrmkMIkkcXc5qFO8Q\nJAmErek750rM7F5gAZABPOOcW2tmE4Dlzrm5ZnYB8ArQHBhhZo8753oA3YDfmdlJAm8wkyqM+hGp\nEZ/JWG2bZsYgkjgr9wlgIf9u8rr92kdDFu9c8KrX8sySuryWYXDOzQfmVzj2SLnvlxFo9ql43dtA\nr1rGKFJtHZo3jHcIcWFU3WmrDl3RjFxJOj4duDltm8QgksTjswWkhm+mNyV9STnpuLBadai2n96U\n9CWp+NTyU21hteryGaaq2n76UtIXSUOq7acvJX1JGj47Y6XLZKxwNHxTKqOkL0nDZ2csCVj4o8Fh\ny2jz9PSkpC8po1EDbQRXXpMGGVWe1+bp6UlJX5KCTwfuhee0jEEkyeODx4eHLZOlDt20o6QvKUEv\nZBE/+luRhOdTy79EHbgh+ay++XDh6hhEIolCSV8kzf1p6bZ4hyAxpKQvCU3DNGtPwzelPCV9SWga\npll7PsM3C1d5b5EhSU5JX5JahmmdHR91w/yaxs4uik0gEndK+pKwXvPowB2Sk97r7Pja/PPwHbo5\n4+eHLSPJT0lfEla49WG0mmZkHTmhFXnSgZK+JKTzHgo/aSjdV9OsLm2eLqCkLwmqRJXOuNDm6alP\nSV8Sjs9IkhaN6scgktQz4NwW8Q5B4kxJXxKOz0iSPp2axyCS1DPrjgvDltF6PKlNSV+SToO6etnW\nRmaGOsDTmf56JKH4dOBe3KV1DCJJXRsmXhG2TO9H/x6DSCQelPQloagDNzbCTdb66uiJ2AQiMaek\nLwnDp3apdXYiw2eyllbfTE3aakgShmqXETRiROjj8+Z53+JPS7fxxKheEQpIEoVq+pI0NEwzssb0\n7xTvECQOlPQlIfgME9QwzcjyqcVrhm7qUdKXpKBNz6Mj3OBNzdBNPfpLkrhbfF5fZoQpo03Po+Pj\nSVdqMlaaUU1f4u5kvANIc+GGb+pNIbUo6UtcDZvyetgyGqYZXT7DNyV1KOlLXG3a9XW8QxAPmqGb\nOpT0JaFpnZ3EoDkUqUN/URI3Pm3FWmcnNnw2WNF2iqnBK+mb2XAz22hmm82sIMT5gWa20sxKzOza\nCuduNrNNwa+bIxW4pD5NxoqtcGvtazvF1BB2yKaZZQBTgWFAMbDMzOY659aVK7YNuAV4sMK1LYBH\ngXwCW56uCF67NzLhS7Ly6cDVZKwoqGx5BmDWvHkaqZMGfGr6fYHNzrmPnHPHgBeBkeULOOe2Ouc+\n4PTRd5cDC51zXwYT/UJgeATiliSnDtzkpBm6yc8n6bcHtpd7XBw85qM210qK0jDNxNXlrEZVntcM\n3eSXEB25ZnanmS03s+W7d++OdzgSZarlJ66FPxoctow6dJObT9LfAXQs97hD8JgPr2udc9Odc/nO\nufzWrTVaQySRqUM3ufkk/WVAFzPrbGb1gdHAXM/7LwAuM7PmZtYcuCx4TNKUT0ehmnbiy2f4ptr2\nk1fYpO+cKwHuJZCs1wN/ds6tNbMJZnY1gJldYGbFwHXA78xsbfDaL4GfEnjjWAZMCB4TCSkh2htF\nq2+mMK9VNp1z84H5FY49Uu77ZQSabkJd+wzwTC1ilBTh0xZ8iWr5CUGrb6YuVawkZtQWnFr0ppCc\nlPQlJvpNXBi2TM92TWMQifjKzAjXyAPnPaTEn2yU9CUmdh44FrZM26aZMYhEfG2YeEXYMiX68JZ0\nlPQl6gpXhR/hq3V2ElObxvp/STXaLlGibuzsIgBmzHm80jJaZycxvTt+WNi2+6yCV72GeUpiUE1f\n4q5D84bxDkGqEG47RUkuqulLVPnsuJTTtkkMIpGwKlmBczOwaP3O047ffu2jZd/3m7iQd8cPi1Zk\nEkGq6UtUhdtxqV4dVSOTQaMGVdcPfTrqJTEo6UvU+NTyB2WfFYNIpLYuPKdl2DI+w3Il/tS8IzVX\n2YYc8+YB2lc13ai2nxxU05eo8Jm0o8lYycWnw91neK7El5K+RIXPpB1NxkouPh3uPwwOz5XEpaQv\ncdGgrl56yShcbd8BDxeujk0wUiP6y5OI81mI6+Iu2iwnGfnU9v+0dFsMIpGaUtKXiPp8/5GwZTQZ\nK7llWPhhtmrbT1xK+hJRaz7dH7aMJmMltyE54YfZqm0/cSnpS0xpYbXUEG5SnRbfTFxK+hIxoabq\nV6SF1VKDz6Q6rbWfmJT0JSL+tWl32DLhpvJLcgm3EJvW2k9MSvoSEUdLToYt4zOVX5LH5p+HX05Z\nWyomHiV9qTWfETtj+neKQSQSa9pkJfko6Uut+YzYeWJUrxhEIrHms5yyavuJRUlfoq7LWY3iHYJE\nkf5/k4s5l1i9Lfn5+W758uXxDkN8jBjhNWJnaLc2MQhGYi64miqEr80b8LG2VIwqM1vhnMsPV07D\nKaTGNnz+Vdgymn2bHto0rl+2tHKleyEr6ScENe9IjRXvPRy2jGbfpgeftn2fTXUk+pT0pUZ8VlLU\n7Nv00qRBRpXntalOYlDSlxrxWUlRs2/TywePDw9bRiN54k9t+lJtWi9dgJDbZT67cRfHw1w2bMrr\nLPzR4KiEJOGppi/V5lPL14id9OSzJs+mXV/HIBKpjJK+VEvO+Plhy6gtP735JBUtxhY/SvpSLUdO\nhJ/Xobb89HaJx6c8LcYWP0r64s2ndqaVNAX8Xged1akbF0r64qVw1Q6v2plW0hTwex2osh8fSvri\nZazH9nc92zWNQSSSLHw68zWEM/a8kr6ZDTezjWa22cwKQpxvYGazg+ffNbOs4PEsMztsZkXBr99G\nNnyJhWFTXvcq17ZpZnQDkaTjswzHjb9/JwaRSKmwDW9mlgFMBYYBxcAyM5vrnFtXrthtwF7n3Hlm\nNhr4b+D64LktzrncCMctMfR/p/04bBkN0ZRQcto2Cbtcx1tbvoxRNAJ+Nf2+wGbn3EfOuWPAi8DI\nCmVGAs8Gv58DXGpmYTZTk2Tg8/FbbYRSlQyPVNBv4sIYRCLg9/faHthe7nFx8FjIMs65EmA/UNqT\n09nMVpnZG2Z2cS3jlRjybdbxGaIn6WtITvgJWzsPHFMzT4xEe3zdZ0An59weMzsfKDSzHs65U9bk\nNbM7gTsBOnXStnqJwmfmZL06+kAn4fkM4VQzT2z41PR3AB3LPe4QPBayjJnVBZoCe5xzR51zewCc\ncyuALUDXik/gnJvunMt3zuW3bt26+j+FRJxvrctn2r2I71Be1fajzyfpLwO6mFlnM6sPjAbmVigz\nF7g5+P21wGLnnDOz1sGOYMzsHKAL8FFkQpdo8ql1qfNWqmOrxyYqqu1HX9jPXM65EjO7F1gAZADP\nOOfWmtkEYLlzbi7wv8BzZrYZ+JLAGwPAQGCCmR0HTgJ3Oef0v5rgfGZKakcsqbYRI/jDhl2cqLBF\n6+3XPnrK45zx89kw8YpYRpZWvNr0nXPzgfkVjj1S7vsjwHUhrvsL8Jdaxigx1LngVa+ZktoRS2pi\nSM5ZYfdVPnLCafnlKNJoOykzbMrrXglfM2+lNnyaBbX8cvQo6UsZnz+0BnXraOat1JrP8tsaux8d\nSvoC+G9afXEXja6S2vNZfltj96ND6+Cmowrb3G34/Cum7D18WodaRWrWkUjq0LyhlmiIA9X009zn\n+4+E/cODwCQsNetIJPkOBtBKnJGlpJ/m1n32VfhCaBKWRIfvXA8l/shR0k9ji9bv5KQLP15Hk7Ak\nmnwmbYH/WlBSNSX9NBVurHSpNo21yblEn8/rbNOurylcVXEFGKkuJf00tNgz4WdmGO+OHxblaETw\nfp390GMHN6maOY+P97GUn5/vli9fHu8wUla/iQuZ+IfxXmXVrCOxVv4TaGWjyeoabP65X5NQOjGz\nFc65/HDlNGQzlVUYmgkw0bOWr7V1JB58hnGWuEDHrm9fgJxKzTtpxLcdv0HdOlpbR+Iip20Tr522\nQCN6akpJP034JvwWjepr1q3E1ZCcs7wTk8+KsHIqJf004Jvw6+A3PV4k2i7p1sZrRI9DQzmrS236\nKc434YP2upXE8u7SX4d8/Vbs4N2062seLlzNE6N6xSq0pKakn8Kqk/A1UkcSkU/HLsCflm5j2I9v\nDz1zfN68KESWvNS8k6Kq09apkTqSqHLaNvF+fR4/6apV0UlXquknkxBDMIHTajI54+d7bYYC0KhB\nXY3UkYSW07YJzRrWZ82n+70Rar01AAAK7ElEQVTKL1q/U59cq6Cafoo576FXOXLCL+XXq2NceE7L\nKEckUnttm2ZW6xOpavyVU00/RQyb8nq1tphr1KCuEr4klZy2TTh07IR3edX4Q1PST0SVNeNUot/E\nhew8cMy7fIO6dZTwJSn16dSc5wp/ytGSk17lF63fyfiJC7WGVDlK+knutfU72dnDP+HXq2OafCVJ\n7eIurfnXpt3eiX/iH8az+A/jTx+SnKajepT0k9Rr63d6d9aW6tmuqXa/kpRwcZfWvPPRHr4+WuJV\n/iSBWn+DunXSvtKjpJ9kPt9/5LRRDDPmPF7lNXXQxCtJPRee05INn3/lNY6/1NGSk7y2fieXpvHf\ng0bvJJGV2/Z6D1srlWGmhC8pK6dtk2p31joCtf6HC1dHJ6gEp/X0E1GFjtxQtXsfLRrV11o6kjbe\n2LiL4yern886NG946lyVJG3r13r6KWDD51+xY+/harfdQ4gXskiKG5R9Vo0Sf/HewxTvPZw2wzuV\n9BNQbSeWpMuLV6SiQdln1fiT8aL1OzHgwKodjMprH/ngEoSSfrxUaMKpbodUKJpwJRKYvdu2aSZL\nNuziRDWbrx1w5nXfYUkdo1vbJqePdkvSpp/ylPSjLcxEq5q2Q1ak5hyRUw3JCay4WZO/sRMnHWs+\n3V/2icGA9s0bkhPpIONAo3fi4PP9R1i8YReL1u+sdcJv0ag+Q7u1UcIXqcSg7LNo1KB29VtHoO3/\nnIdepXDVjsgEFieq6UfZ5/uPsHn3QY4c918zxEeGWVlNRkSqVpMx/aGcdDB2dhFjZxcBMODcFsy6\n48JIhBgzSvoR9nDhap5/dxulFfgZNehQqkpmvQzOa32mZtaKVFNO2yZln4gjtQrnW1u+JKvgVZo1\nrMdjV/dIig7g9Bmn77kWvY/CVTt4bO5a9h0+Xnbs2Vd+Sh3Dez2Q6tKYe5HIq0lnbyi3X/so9eoY\n1/ftyJINu/l032HaNWvIuMuzY/ZGENFx+mY2HPgVkAHMcM5NqnC+AfBH4HxgD3C9c25r8NxDwG3A\nCeB+59yCavwcMfFw4WpeeHc7J5zDgDoGv3sp9NIGFffnLHX8RHSSfb06FnoLOBGptSE5Z7Fy216+\n/Np/0cLKHD/pmLV0W9m8mh37DvPQy4FZv6Py2lO4ageTF2wM+YZQcaXcNo3rR21l0LBJ38wygKnA\nMKAYWGZmc51z68oVuw3Y65w7z8xGA/8NXG9m3YHRQA+gHbDIzLo65yLbwA1V/kKBSv9jx1f4ZTvA\ncw+SqFETjkjsVPwEXZ2F3CqqmDoOHz/B5AUbAXjo5dUcDvbtPT6jgIxnjM/bNmHTrgNMDLYQlFYq\ndx44Rr8oLQntU9PvC2x2zn0EYGYvAiOB8kl/JPBY8Ps5wNNmZsHjLzrnjgIfm9nm4P3eiUz4AYWr\ndpzyC634Dnvj79/hB5W8k0/8w/hqPVe4xc1qQ004IvFXOtelOss3V+XTfYeZvGBjWX4qdeKkY/Pu\ng5U+R3X2yKgOn6TfHthe7nEx0K+yMs65EjPbD7QMHl9a4dqIN3CV/4WWT8qZczLgvFb8IEG3TtNk\nKpHEVX4J5s/3H2HD519RUskQ66oqg5n1Mjh6/ETI5VSORnhUn4+EGL1jZncCdwJ06tSp2td/uu/f\nw7DKt7kb8PGkK7m94NVaxxgJ0WynE5HoaRv8KlW4agePz1vL3kOBwRzNGtbjqm9+g7+s2HFKjb5h\nvQx+fk0vJi/YyI59pw8Xbd+sYcjj0eST9HcAHcs97hA8FqpMsZnVBZoS6ND1uRbn3HRgOgRG7/gG\nX6pdJb+4ds38N1KurWYN69GjXWOWfrSXE86RYcYN/TryxKheMYtBRGJjVF77kKNy8s9uUWnfYvkm\naAi8IYy7PJufz18XsimnTeP6UYndJ+kvA7qYWWcCCXs08P0KZeYCNxNoq78WWOycc2Y2F3jezKYQ\n6MjtArwXqeBLjbs8u9JfKAQmULy15ctKr2/TuD5fHDx+yuidqjpz28d4KJaIJIfK3gxKj4V6QxiV\n1z6mo3e8xumb2RXAUwSGbD7jnJtoZhOA5c65uWaWCTwH5AFfAqPLdfyOB24FSoCxzrm/VfVcNR2n\nH270zo2/fydk4k/GGXUiIhX5jtNPn8lZIiIpzDfpa8E1EZE0oqQvIpJGlPRFRNKIkr6ISBpR0hcR\nSSNK+iIiaURJX0QkjSjpi4ikESV9EZE0oqQvIpJGlPRFRNKIkr6ISBpR0hcRSSNK+iIiaURJX0Qk\njSjpi4ikESV9EZE0knA7Z5nZbuCTKNy6FfBFFO4bTckWs+KNLsUbXckWL5wa89nOudbhLki4pB8t\nZrbcZyuxRJJsMSve6FK80ZVs8ULNYlbzjohIGlHSFxFJI+mU9KfHO4AaSLaYFW90Kd7oSrZ4oQYx\np02bvoiIpFdNX0Qk7aVl0jezB8zMmVmreMdSFTObbGYbzOwDM3vFzJrFO6ZQzGy4mW00s81mVhDv\neKpiZh3NbImZrTOztWb2X/GOyYeZZZjZKjP7a7xj8WFmzcxsTvD1u97MLox3TFUxsx8GXw9rzOwF\nM8uMd0zlmdkzZrbLzNaUO9bCzBaa2abgv8197pV2Sd/MOgKXAdviHYuHhUBP51xv4EPgoTjHcxoz\nywCmAt8GugM3mFn3+EZVpRLgAedcd6A/8J8JHm+p/wLWxzuIavgV8HfnXA7wTRI4djNrD9wP5Dvn\negIZwOj4RnWamcDwCscKgNecc12A14KPw0q7pA88CfwYSPjODOfcP5xzJcGHS4EO8YynEn2Bzc65\nj5xzx4AXgZFxjqlSzrnPnHMrg98fIJCM2sc3qqqZWQfgSmBGvGPxYWZNgYHA/wI454455/bFN6qw\n6gINzawucAbwaZzjOYVz7p/AlxUOjwSeDX7/LDDK515plfTNbCSwwzn3frxjqYFbgb/FO4gQ2gPb\nyz0uJsGTaCkzywLygHfjG0lYTxGoqJyMdyCeOgO7gT8Em6RmmFmjeAdVGefcDuCXBD79fwbsd879\nI75ReWnjnPss+P3nQBufi1Iu6ZvZomC7XMWvkcBPgEfiHWN5YeItLTOeQLPErPhFmlrM7EzgL8BY\n59xX8Y6nMmZ2FbDLObci3rFUQ12gD/Ab51we8DWeTQ/xEGwLH0ngzaod0MjMxsQ3qupxgWGYXq0X\ndaMcS8w554aGOm5mvQj8p75vZhBoKllpZn2dc5/HMMRTVBZvKTO7BbgKuNQl5vjaHUDHco87BI8l\nLDOrRyDhz3LOvRzveMIYAFxtZlcAmUATM/uTcy6Rk1IxUOycK/0ENYcETvrAUOBj59xuADN7GfgW\n8Ke4RhXeTjP7hnPuMzP7BrDL56KUq+lXxjm32jl3lnMuyzmXReCF2SeeCT8cMxtO4GP91c65Q/GO\npxLLgC5m1tnM6hPoAJsb55gqZYF3/P8F1jvnpsQ7nnCccw855zoEX7OjgcUJnvAJ/k1tN7Ps4KFL\ngXVxDCmcbUB/Mzsj+Pq4lATueC5nLnBz8Pubgf/nc1HK1fRTzNNAA2Bh8NPJUufcXfEN6VTOuRIz\nuxdYQGDUwzPOubVxDqsqA4D/A6w2s6LgsZ845+bHMaZUdB8wK1gR+Aj4QZzjqZRz7l0zmwOsJNCM\nuooEm51rZi8Ag4FWZlYMPApMAv5sZrcRWJn4e173SswWAxERiYa0ad4RERElfRGRtKKkLyKSRpT0\nRUTSiJK+iEgaUdIXEUkjSvoiImlESV9EJI38f7FSssZOiQOOAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RRUK6UCjKa4Y",
        "colab_type": "text"
      },
      "source": [
        "#### 二维Gibbs采样实例python实现  \n",
        "\n",
        "假设我们要采样的是一个二维正态分布 $N(\\mu, \\Sigma)$ ，其中： $\\mu=(\\mu_{1}, \\mu_{2})= (5, -1)$ , $\\Sigma = \\begin{pmatrix}\n",
        "\\sigma^{2}_{1} &   \\rho \\sigma_{1}\\sigma_{2}b\\rho \\sigma_{2}& \n",
        "\\sigma^{2}_{2}\\end{pmatrix} = \\begin{pmatrix}\n",
        " 1& 1b1 & \n",
        "4\\end{pmatrix}$;\n",
        "\n",
        "而采样过程中的需要的状态转移条件分布为：\n",
        "\n",
        "$P(x_{1}|x_{2}) = N(\\mu_{1}+ \\rho \\sigma_{1}/\\sigma_{2}(x_{2} - \\mu_{2}), (1 - \\rho^{2})\\sigma^{2}_{1})$\n",
        "\n",
        "$P(x_{2}|x_{1}) = N(\\mu_{2}+ \\rho \\sigma_{2}/\\sigma_{1}(x_{1} - \\mu_{1}), (1 - \\rho^{2})\\sigma^{2}_{2})$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3ZSOvaDcMHfv",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 335
        },
        "outputId": "0c3d6e18-ba4e-4e06-b61f-2585cc78715d"
      },
      "source": [
        "from mpl_toolkits.mplot3d import Axes3D\n",
        "from scipy.stats import multivariate_normal\n",
        "\n",
        "samplesource = multivariate_normal(mean=[5,-1], cov=[[1,0.5],[0.5,2]])\n",
        "\n",
        "def p_ygivenx(x, m1, m2, s1, s2):\n",
        "    return (random.normalvariate(m2 + rho * s2 / s1 * (x - m1), math.sqrt(1 - rho ** 2) * s2))\n",
        "\n",
        "def p_xgiveny(y, m1, m2, s1, s2):\n",
        "    return (random.normalvariate(m1 + rho * s1 / s2 * (y - m2), math.sqrt(1 - rho ** 2) * s1))\n",
        "\n",
        "N = 5000\n",
        "K = 20\n",
        "x_res = []\n",
        "y_res = []\n",
        "z_res = []\n",
        "m1 = 5\n",
        "m2 = -1\n",
        "s1 = 1\n",
        "s2 = 2\n",
        "\n",
        "rho = 0.5\n",
        "y = m2\n",
        "\n",
        "for i in range(N):\n",
        "    for j in range(K):\n",
        "        x = p_xgiveny(y, m1, m2, s1, s2)   #y给定得到x的采样\n",
        "        y = p_ygivenx(x, m1, m2, s1, s2)   #x给定得到y的采样\n",
        "        z = samplesource.pdf([x,y])\n",
        "        x_res.append(x)\n",
        "        y_res.append(y)\n",
        "        z_res.append(z)\n",
        "\n",
        "num_bins = 50\n",
        "plt.hist(x_res, num_bins, normed=1, facecolor='green', alpha=0.5,label='x')\n",
        "plt.hist(y_res, num_bins, normed=1, facecolor='red', alpha=0.5,label='y')\n",
        "plt.title('Histogram')\n",
        "plt.legend()\n",
        "plt.show()"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.6/dist-packages/matplotlib/axes/_axes.py:6521: MatplotlibDeprecationWarning: \n",
            "The 'normed' kwarg was deprecated in Matplotlib 2.1 and will be removed in 3.1. Use 'density' instead.\n",
            "  alternative=\"'density'\", removal=\"3.1\")\n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEICAYAAACzliQjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAG+JJREFUeJzt3X+UHWWd5/H3xyYhu8DwI+lRTCek\nkWbWgByiTVRmxAiBBGdMYNZZGsgBFjwRJTPOeNwRFja2YZiDuMsBj1HIaA7G0Qk/XN1eNxhRiHNc\nCXYjPxPIpgmBdIYfmQ6/IiSQ5Lt/VDVTudzuru6+9/btrs/rnD5d9dTz1P3e6tvf+9yn6j6liMDM\nzIrhXaMdgJmZ1Y6TvplZgTjpm5kViJO+mVmBOOmbmRWIk76ZWYE46du4JGmDpDmjHYdZvXHStzFJ\n0lZJc0vKLpH0a4CIOCEi1g2yjxmSQtJBVQzVrK446ZtVid9MrB456du4lP0kIGm2pC5Jr0p6QdKN\nabV/Tn+/LGmXpI9KepekayQ9I+lFSaskHZ7Z70Xptl5J/63kcdol3SXpHyW9ClySPvb9kl6W9Jyk\nb0qamNlfSPq8pM2SXpN0raT3SfpNGu8d2fpmI+Wkb0VwM3BzRPwB8D7gjrT8tPT3ERFxaETcD1yS\n/nwCOBY4FPgmgKSZwLeAC4GjgcOBqSWPtRC4CzgC+AGwD/gbYArwUeAM4PMlbeYBHwI+AvwtsAJY\nBEwDTgTOH8FzNzuAk76NZT9Je9AvS3qZJCGX8xZwnKQpEbErItYPsM8LgRsjYktE7AKuAtrSoZpP\nA/87In4dEW8CS4HSyavuj4ifRMT+iHgjIh6MiPURsTcitgK3Ah8vaXNDRLwaERuAx4Gfp4//CnA3\nMCv/ITEbmJO+jWXnRMQRfT+8swfd5zLgeOBJSZ2S/myAfb4XeCaz/gxwEPDudNu2vg0R8TrQW9J+\nW3ZF0vGSfirp+XTI5+9Jev1ZL2SW3yizfugA8ZoNiZO+jXsRsTkizgf+EPgacJekQ3hnLx3gX4Bj\nMuvTgb0kifg5oKlvg6R/B0wufbiS9W8DTwIt6fDSfwU0/GdjNjJO+jbuSVokqTEi9gMvp8X7gR3p\n72Mz1f8J+BtJzZIOJemZ3x4Re0nG6j8l6dT05Go7gyfww4BXgV2S/gPwuUo9L7PhcNK3IpgPbJC0\ni+Skbls63v46cB3wf9PzAh8BVgLfJ7my52lgN/CXAOmY+18Cq0l6/buAF4E9Azz2l4ALgNeAfwBu\nr/zTM8tPvomK2fCknwReJhm6eXq04zHLwz19syGQ9ClJ/z49J/DfgceAraMblVl+TvpmQ7OQ5GTv\nvwAtJENF/rhsY4aHd8zMCsQ9fTOzAqm7CaGmTJkSM2bMGO0wzMzGlAcffPBfI6JxsHp1l/RnzJhB\nV1fXaIdhZjamSHpm8Foe3jEzKxQnfTOzAnHSNzMrkLob0zczGy1vvfUWPT097N69e7RD6dekSZNo\nampiwoQJw2rvpG9mlurp6eGwww5jxowZSPU3GWpE0NvbS09PD83NzcPah4d3zMxSu3fvZvLkyXWZ\n8AEkMXny5BF9EnHSNzPLqNeE32ek8eVK+pLmS9okqVvSlQPU+4/pjZ5bM2VXpe02SZo3omjNzGxE\nBh3Tl9QALAfOBHqATkkdEbGxpN5hwBeABzJlM4E24ASSW839QtLxEbGvck/BzKw62te1V3Z/cyq7\nv+HIcyJ3NtAdEVsAJK0mmWlwY0m9a0luRfdfMmULgdURsQd4WlJ3ur/7Rxq4mY1N2URaD0mwaPIM\n70zlwJs996Rlb5P0QWBaRPyfobZN2y+W1CWpa8eOHbkCN7Oxo31d+9s/NrDOzk5OOukkdu/eze9/\n/3tOOOEEHn/88Yrtf8SXbEp6F3AjcMlw9xERK4AVAK2trZ7r2cwK65RTTmHBggVcc801vPHGGyxa\ntIgTTzyxYvvPk/S3A9My601pWZ/DgBOBdelZ5fcAHZIW5GhrZmYlli5dyimnnMKkSZP4xje+UdF9\n5xne6QRaJDVLmkhyYrajb2NEvBIRUyJiRkTMANYDCyKiK63XJulgSc0kdxr6bUWfgZnZONPb28uu\nXbt47bXXKv7t4EGTfkTsBZYAa4EngDsiYoOkZWlvfqC2G4A7SE76/gy4wlfumJkN7LOf/SzXXnst\nF154IV/+8pcruu9cY/oRsQZYU1K2tJ+6c0rWrwOuG2Z8ZmajZjSuLlq1ahUTJkzgggsuYN++fZx6\n6qnce++9nH766RXZv+feMTOrIxdddBEXXXQRAA0NDTzwwAODtBgaJ30zqwpfnlmfPPeOmVmBOOmb\nmRWIk76ZWYE46ZuZFYiTvplZgfjqHTOrmKFesdNXv25n22xvr+/9DYN7+mZmBeKkb2ZWR5YuXcpN\nN9309vrVV1/NzTffXLH9O+mbmdWRSy+9lFWrVgGwf/9+Vq9ezaJFiyq2f4/pm5nVkRkzZjB58mQe\neughXnjhBWbNmsXkyZMrtn8nfTOzOvOZz3yG2267jeeff55LL720ovv28I6ZWZ0599xz+dnPfkZn\nZyfz5s2r6L7d0zcz688oXWI5ceJEPvGJT3DEEUfQ0NBQ0X076ZuZ1Zn9+/ezfv167rzzzorvO9fw\njqT5kjZJ6pZ0ZZntl0t6TNLDkn4taWZaPkPSG2n5w5JuqfQTMDMbTzZu3Mhxxx3HGWecQUtLS8X3\nP2hPX1IDsBw4E+gBOiV1RMTGTLUfRsQtaf0FwI3A/HTbUxFxcmXDNjMbn2bOnMmWLVuqtv88Pf3Z\nQHdEbImIN4HVwMJshYh4NbN6CBCVC9HMrHYi6jt9jTS+PEl/KrAts96Tlh1A0hWSngJuAP4qs6lZ\n0kOSfiXpY+UeQNJiSV2Sunbs2DGE8M3MKmfSpEn09vbWbeKPCHp7e5k0adKw91GxE7kRsRxYLukC\n4BrgYuA5YHpE9Er6EPATSSeUfDIgIlYAKwBaW1vr82ib2bjX1NRET08P9dz5nDRpEk1NTcNunyfp\nbwemZdab0rL+rAa+DRARe4A96fKD6SeB44GuYUVrZlZFEyZMoLm5ebTDqKo8Sb8TaJHUTJLs24AL\nshUktUTE5nT1T4HNaXkjsDMi9kk6FmgBqneGwsxqrhI3QM/uo26nWR4nBk36EbFX0hJgLdAArIyI\nDZKWAV0R0QEskTQXeAt4iWRoB+A0YJmkt4D9wOURsbMaT8TMzAaXa0w/ItYAa0rKlmaWv9BPux8B\nPxpJgGZmVjmee8fMrECc9M3MCsRJ38ysQJz0zcwKxEnfzKxAnPTNzArESd/MrECc9M3MCsRJ38ys\nQJz0zcwKxEnfzKxAnPTNzArESd/MrECc9M3MCsRJ38ysQJz0zcwKpGI3Rjez4qjELRLz7Nu3Tqy8\nXD19SfMlbZLULenKMtsvl/SYpIcl/VrSzMy2q9J2myTNq2TwZmY2NIMmfUkNwHLgbGAmcH42qad+\nGBEfiIiTgRuAG9O2M0lupH4CMB/4Vro/MzMbBXl6+rOB7ojYEhFvAquBhdkKEfFqZvUQINLlhcDq\niNgTEU8D3en+zMxsFOQZ058KbMus9wAfLq0k6Qrgi8BE4PRM2/UlbaeWabsYWAwwffr0PHGbmdkw\nVOzqnYhYHhHvA74MXDPEtisiojUiWhsbGysVkpmZlciT9LcD0zLrTWlZf1YD5wyzrZmZVVGepN8J\ntEhqljSR5MRsR7aCpJbM6p8Cm9PlDqBN0sGSmoEW4LcjD9vMzIZj0DH9iNgraQmwFmgAVkbEBknL\ngK6I6ACWSJoLvAW8BFyctt0g6Q5gI7AXuCIi9lXpuZiZ2SByfTkrItYAa0rKlmaWvzBA2+uA64Yb\noJmZVY6nYTAzKxAnfTOzAnHSNzMrECd9M7MCcdI3MysQJ30zswJx0jczKxAnfTOzAnHSNzMrECd9\nM7MCcdI3MysQJ30zswLJNeGamdloaF/X/m/Lc9r7rWf5uadvZlYgTvpmZgXipG9mViAe0zez3LJj\n7DY25erpS5ovaZOkbklXltn+RUkbJT0q6ZeSjsls2yfp4fSno7StmZnVzqA9fUkNwHLgTKAH6JTU\nEREbM9UeAloj4nVJnwNuAM5Lt70RESdXOG4zMxuGPD392UB3RGyJiDeB1cDCbIWIuC8iXk9X1wNN\nlQ3TzMwqIU/Snwpsy6z3pGX9uQy4O7M+SVKXpPWSzinXQNLitE7Xjh07coRkZmbDUdETuZIWAa3A\nxzPFx0TEdknHAvdKeiwinsq2i4gVwAqA1tbWqGRMZmb2b/L09LcD0zLrTWnZASTNBa4GFkTEnr7y\niNie/t4CrANmjSBeMzMbgTxJvxNokdQsaSLQBhxwFY6kWcCtJAn/xUz5kZIOTpenAH8MZE8Am5lZ\nDQ06vBMReyUtAdYCDcDKiNggaRnQFREdwNeBQ4E7JQE8GxELgPcDt0raT/IGc33JVT9mZlZDucb0\nI2INsKakbGlmeW4/7X4DfGAkAZqZWeV4GgYzswJx0jczKxAnfTOzAnHSNzMrECd9M7MC8dTKZjYm\n+NaJleGevplZgTjpm5kViJO+mVmBOOmbmRWIk76ZWYE46ZuZFYgv2TSzAWUvlbSxzz19M7MCcdI3\nMysQJ30zswJx0jczK5BcSV/SfEmbJHVLurLM9i9K2ijpUUm/lHRMZtvFkjanPxdXMngzMxuaQZO+\npAZgOXA2MBM4X9LMkmoPAa0RcRJwF3BD2vYo4CvAh4HZwFckHVm58M3MbCjy9PRnA90RsSUi3gRW\nAwuzFSLivoh4PV1dDzSly/OAeyJiZ0S8BNwDzK9M6GZmNlR5rtOfCmzLrPeQ9Nz7cxlw9wBtp5Y2\nkLQYWAwwffr0HCGZDVF7e74ys3GuoidyJS0CWoGvD6VdRKyIiNaIaG1sbKxkSGZmlpGnp78dmJZZ\nb0rLDiBpLnA18PGI2JNpO6ek7brhBGqWi3vvZgPK09PvBFokNUuaCLQBHdkKkmYBtwILIuLFzKa1\nwFmSjkxP4J6VlpmZ2SgYtKcfEXslLSFJ1g3AyojYIGkZ0BURHSTDOYcCd0oCeDYiFkTETknXkrxx\nACyLiJ1VeSZmQ9XfpwJ/Wqh7vnXi8OWacC0i1gBrSsqWZpbnDtB2JbByuAGa1ZxP+to45lk2bexy\nIjYbMk/DYGZWIE76ZmYF4qRvZlYgTvpmZgXipG9mViBO+mZmBeJLNq3++dJMs4px0jfLo2Df3s1+\n49XGFw/vmJkViJO+mVmBOOmbmRWIk76ZWYE46ZuZFYiTvplZgTjpm5kVSK7r9CXNB24muXPWdyLi\n+pLtpwE3AScBbRFxV2bbPuCxdPXZiFhQicBtHBqL17wX7Pp9G/sGTfqSGoDlwJlAD9ApqSMiNmaq\nPQtcAnypzC7eiIiTKxCrmZmNUJ6e/mygOyK2AEhaDSwE3k76EbE13ba/CjGamVmF5BnTnwpsy6z3\npGV5TZLUJWm9pHOGFJ2ZmVVULebeOSYitks6FrhX0mMR8VS2gqTFwGKA6dOn1yAkM7NiytPT3w5M\ny6w3pWW5RMT29PcWYB0wq0ydFRHRGhGtjY2NeXdtZmZDlKen3wm0SGomSfZtwAV5di7pSOD1iNgj\naQrwx8ANww3WzKxUdkbQ9jnt/dazxKA9/YjYCywB1gJPAHdExAZJyyQtAJB0iqQe4C+AWyVtSJu/\nH+iS9AhwH3B9yVU/ZmZWQ7nG9CNiDbCmpGxpZrmTZNintN1vgA+MMEYzM6sQfyPXzKxAfOcsGx3+\nxmrd8d2yisE9fTOzAnFP36wayn2S8acbqwPu6ZuZFYiTvplZgTjpm5kViJO+mVmBOOmbmRWIk76Z\nWYE46ZuZFYiTvplZgTjpm5kViL+Ra1Yr/X0j19/UrZi++YM8r37/nPStupzQzOqKh3fMzArESd/M\nrEByDe9Img/cDDQA34mI60u2nwbcBJwEtEXEXZltFwPXpKt/FxHfq0TgZjZynkO/eAbt6UtqAJYD\nZwMzgfMlzSyp9ixwCfDDkrZHAV8BPgzMBr6S3izdzMxGQZ7hndlAd0RsiYg3gdXAwmyFiNgaEY8C\n+0vazgPuiYidEfEScA8wvwJxm5nZMORJ+lOBbZn1nrQsj1xtJS2W1CWpa8eOHTl3bWZmQ1UXJ3Ij\nYkVEtEZEa2Nj42iHY2Y2buVJ+tuBaZn1prQsj5G0NTOzCstz9U4n0CKpmSRhtwEX5Nz/WuDvMydv\nzwKuGnKUNjb4i1hmdW/QpB8ReyUtIUngDcDKiNggaRnQFREdkk4BfgwcCXxK0lcj4oSI2CnpWpI3\nDoBlEbGzSs/FbGzyTdSthnJdpx8Ra4A1JWVLM8udJEM35dquBFaOIEYzM6uQujiRa2ZmteEJ18xs\n3Ml+09gzbh7IPX0zswJx0jczKxAnfTOzAnHSNzMrECd9M7MC8dU7ZvWoivfT9Rz6xeakb0Pnb4ua\njVke3jEzKxAnfTOzAnHSNzMrECd9M7MCcdI3MysQJ30zswLxJZtmNq55xs0DuadvZlYguXr6kuYD\nN5PcLvE7EXF9yfaDgVXAh4Be4LyI2CppBvAEsCmtuj4iLq9M6FZ1/hJW/aniN3WtGAZN+pIagOXA\nmUAP0CmpIyI2ZqpdBrwUEcdJagO+BpyXbnsqIk6ucNxmNkSefsEg3/DObKA7IrZExJvAamBhSZ2F\nwPfS5buAMySpcmGamVkl5En6U4FtmfWetKxsnYjYC7wCTE63NUt6SNKvJH2s3ANIWiypS1LXjh07\nhvQEzMwsv2qfyH0OmB4Rs4AvAj+U9AellSJiRUS0RkRrY2NjlUMyMyuuPEl/OzAts96UlpWtI+kg\n4HCgNyL2REQvQEQ8CDwFHD/SoM3MbHjyXL3TCbRIaiZJ7m3ABSV1OoCLgfuBTwP3RkRIagR2RsQ+\nSccCLcCWikVvleOrP8wKYdCkHxF7JS0B1pJcsrkyIjZIWgZ0RUQH8F3g+5K6gZ0kbwwApwHLJL0F\n7Acuj4id1XgiZmY2OEXEaMdwgNbW1ujq6hrtMIrHPf3xp+Rv6ks2DzTevp0r6cGIaB2snr+Ra2ZW\nIJ57x2wcc+/eSrmnb2ZWIE76ZmYF4uGdovEJ2+Job2fO1nXvKF53yZyah2L1w0nfzAqpqPPse3jH\nzKxA3NM3G2fWlRnSMevjnr6ZWYG4pz+e+aStlTHntnXvKPPJ3eJw0jezwivSSV0nfbNxwOP4lpeT\n/njgYRwboXJDPuBhn/HIJ3LNzArEPX2zMcpDOjYcTvpjjYdyrIaKOOwz3k/qOumb2ZAV5bLPvjeA\n8ZT8cyV9SfOBm0lul/idiLi+ZPvBwCrgQ0AvcF5EbE23XQVcBuwD/ioi1lYs+vHMPXorw0M6NlKD\n3i5RUgPw/4AzgR6SG6WfHxEbM3U+D5wUEZdLagPOjYjzJM0E/gmYDbwX+AVwfETs6+/xxvXtEp3I\nbRjGY6Ify58K6rXXn/d2iXl6+rOB7ojYku54NbAQ2JipsxBoT5fvAr4pSWn56ojYAzyd3jh9NnB/\n3idS95zIrYLGY4Ivp79zBeXU2xtEnruR1esbA+RL+lOBbZn1HuDD/dWJiL2SXgEmp+XrS9pOLX0A\nSYuBxenqLkmbcsQ1BfjXHPVGQz3HBo5vJOo5NhiP8X3vV9WJ5J0qduy+ylcrsZtSg8V3TJ6d1MWJ\n3IhYAawYShtJXXk+yoyGeo4NHN9I1HNs4PhGop5jg8rFl+fLWduBaZn1prSsbB1JBwGHk5zQzdPW\nzMxqJE/S7wRaJDVLmgi0AR0ldTqAi9PlTwP3RnKGuANok3SwpGagBfhtZUI3M7OhGnR4Jx2jXwKs\nJblkc2VEbJC0DOiKiA7gu8D30xO1O0neGEjr3UFy0ncvcMVAV+4M0ZCGg2qsnmMDxzcS9RwbOL6R\nqOfYoELxDXrJppmZjR+ecM3MrECc9M3MCqRuk76kv5C0QdJ+Sa0l266S1C1pk6R5/bRvlvRAWu/2\n9CR0tWK9XdLD6c9WSQ/3U2+rpMfSejX72rGkdknbMzF+sp9689Nj2i3pyhrG93VJT0p6VNKPJR3R\nT72aHb/BjkV6ccLt6fYHJM2oZjwljz1N0n2SNqb/I18oU2eOpFcyf/OltYovffwB/1ZKfCM9fo9K\n+mCN4vqjzDF5WNKrkv66pE5Nj52klZJelPR4puwoSfdI2pz+PrKfthendTZLurhcnXeIiLr8Ad4P\n/BGwDmjNlM8EHgEOBpqBp4CGMu3vANrS5VuAz9Uo7v8BLO1n21Zgyigcy3bgS4PUaUiP5bHAxPQY\nz6xRfGcBB6XLXwO+NprHL8+xAD4P3JIutwG31/DveTTwwXT5MJJpUkrjmwP8tNavtbx/K+CTwN2A\ngI8AD4xCjA3A88Axo3nsgNOADwKPZ8puAK5Ml68s9z8BHAVsSX8fmS4fOdjj1W1PPyKeiIhy38x9\ne2qHiHga6Jva4W3pFBCnk0wJAfA94Jxqxpt53P9EMt/QWPP2dBsR8SbQN91G1UXEzyNib7q6nuT7\nHKMpz7FYSPK6guR1dkb696+6iHguIn6XLr8GPEGZb7rXuYXAqkisB46QdHSNYzgDeCoinqnx4x4g\nIv6Z5KrHrOzrq7/8NQ+4JyJ2RsRLwD3A/MEer26T/gDKTQtR+oKfDLycSSRlp3+ogo8BL0TE5n62\nB/BzSQ+mU0/U0pL0Y/TKfj4q5jmutXApSQ+wnFodvzzH4oCpR4C+qUdqKh1WmgU8UGbzRyU9Iulu\nSSfUNLDB/1b18Hpro/8O2mgeO4B3R8Rz6fLzwLvL1BnWMRzVaRgk/QJ4T5lNV0fE/6p1PAPJGev5\nDNzL/5OI2C7pD4F7JD2ZvstXNT7g28C1JP+I15IMQV1aicfNK8/xk3Q1yfc5ftDPbqp2/MYiSYcC\nPwL+OiJeLdn8O5Jhi13pOZyfkHw5slbq+m+VnuNbAFxVZvNoH7sDRERIqti19aOa9CNi7jCa5Zna\noZfk4+JBaS9sxNM/DBarkukn/pzkngL97WN7+vtFST8mGUaoyD9C3mMp6R+An5bZVNUpM3Icv0uA\nPwPOiHTAssw+qnb8Sgxl6pEeHTj1SE1ImkCS8H8QEf+zdHv2TSAi1kj6lqQpEVGTydhy/K1Ge4qW\ns4HfRcQLpRtG+9ilXpB0dEQ8lw57vVimznaS8w99mkjOgQ5oLA7vDDq1Q5o07iOZEgKSKSKq/clh\nLvBkRPSU2yjpEEmH9S2TnLx8vFzdSisZKz23n8fNM91GteKbD/wtsCAiXu+nTi2P30imHqm69NzB\nd4EnIuLGfuq8p+8cg6TZJP/rNXlTyvm36gAuSq/i+QjwSmY4oxb6/VQ+mscuI/v66i9/rQXOknRk\nOmR7Vlo2sFqdoR7GGe1zScao9gAvAGsz264mubpiE3B2pnwN8N50+ViSN4Nu4E7g4CrHextweUnZ\ne4E1mXgeSX82kAxr1OpYfh94DHg0fTEdXRpfuv5JkitBnqpxfN0kY5MPpz+3lMZX6+NX7lgAy0je\nmAAmpa+r7vR1dmwNj9efkAzVPZo5Zp8ELu97DQJL0uP0CMnJ8VNrGF/Zv1VJfAKWp8f3MTJX6NUg\nvkNIkvjhmbJRO3Ykbz7PAW+lOe8ykvNDvwQ2k9x86qi0bivJ3Qv72l6avga7gf+c5/E8DYOZWYGM\nxeEdMzMbJid9M7MCcdI3MysQJ30zswJx0jczKxAnfTOzAnHSNzMrkP8Pp1vkThLpx98AAAAASUVO\nRK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qmg5nG9J7reE",
        "colab_type": "text"
      },
      "source": [
        "### 习题19.8 "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VEu6ZHpX7wEs",
        "colab_type": "text"
      },
      "source": [
        ""
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OfldxGWecECT",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}